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1.  INTRODUCTION 


The  ballistic  test  ranges  at  the  Arnold  Engineering  Development  Center  (AEDC)  use  two-stage  light- 
gas  guns  to  launch  projectiles  up  to  speeds  of  6  km/s.  Graphite-tipped  models,  with  2.5-in.  diameters 
and  weighing  near  500  g,  are  subjected  to  accelerations  up  to  100,000  g,  pressures  over  40,000  psi, 
and  strain  rates  over  1000  s*1.  More  massive  models  are  subjected  to  even  higher  pressures.  Under 
these  high  loads  and  strain  rates,  die  dynamical  stress  can  be  quite  high.  Occasionally  a  model  will 
fail,  causing  damage  to  the  track  and  scattering  model  fragments  on  die  track.  A  model  failure  not 
only  results  in  die  loss  of  test  data,  it  also  leads  to  expensive  downtime  to  inspect  and  correct  the 
damage. 

The  projectiles  are  often  configured  using  a  model  and  a  sabot  to  hold  the  model  in  place.  They  are 
constructed  from  metals,  plastics,  and  composites.  Under  these  extreme  loads,  yielding  can  occur  with 
the  material  flowing  plastically  at  very  high  strain  rates. 

Although  the  major  loads  act  along  the  axis,  side  loads  are  also  important  to  the  performance  of  the 
gun.  The  model/sabot  must  act  as  a  seal  to  prevent  "blow-by,"  the  high-pressure  gas  leaking  to  the 
front  side  of  the  model.  The  relatively  snug  fit  of  the  projectile  in  die  bore  also  prevents  the  model 
from  ballating  as  it  travels  down  the  barrel.  The  friction  between  the  projectile  and  the  gun-barrel  wall 
tends  to  wear  the  model,  reducing  the  effectiveness  of  the  seal,  increasing  the  side  forces  due  to  gas 
pressure,  and  allowing  more  freedom  for  the  unwanted  ballating  motion. 

The  design  of  the  model/sabot  and  the  gas  gun  attempts  to  satisfy  dual  and  often  opposing 
requirements:  maximizing  the  launch  velocity  while  keeping  the  projectile  intact.  An  understanding  of 
the  dynamics  of  the  projectile  and  its  interaction  with  the  light  gas  and  the  barrel  wall  could  lead  to 
improved  model/sabot  designs  and  procedures  to  reduce  the  possibility  of  model  failure  without 
degrading  the  overall  performance  of  the  gas  gun.  However,  due  to  the  extremely  high  loadings  and 
short  operation  times,  it  is  difficult  at  best  to  measure  the  dynamical  stresses  acting  within  the 
projectile.  Numerical  simulations  of  the  dynamical  process  have  the  potential  to  provide  more  insight 
into  the  complicated  problem  and  offer  a  means  to  judge  strategies  for  improving  the  model 
survivability  and  increasing  the  gas  gun  performance.  It  is  the  development  of  such  a  numerical 
dynamical  stress  model  that  is  the  subject  of  this  effort. 

In  the  discussion  to  follow,  several  key  factors  are  identified  that  are  important  to  the  accuracy  and 
utPty  of  a  numerical  simulation  of  the  dynamics  of  the  model/sabot.  These  include: 

•  The  coupled  response  of  the  model/sabot  and  the  light  gas. 

•  A  good  representation  of  the  particular  geometry  and  construction  of  die  model/sabot. 

•  A  viscoelastic-viscoplastic,  rate-dependent  material  model  that  accounts  for  finite  strains  and 
large  deformations. 

•  A  side-friction  and  normal-force  model  resulting  from  contact  with  the  barrel  wall  or  high- 
pressure  gas  leaking  to  the  sides. 
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During  this  Phase  I  study  conducted  by  QUEST  Integrated,  Inc.,  the  objectives  of  the  effort  have  been 
to  examine  and  develop  models  for  each  of  these  key  areas  and  demonstrate  their  relative  importance. 
In  this  final  report  for  Phase  I,  we  present  the  various  parts  of  the  mathematical  model  and  the  results 
of  some  example  simulations.  The  appendices  contain  the  User's  Guide  for  the  FEAMOD  code,  which 
is  essential  to  this  model,  and  some  file  format  information  and  example  files.  First,  however,  in 
Section  2  we  give  some  background  information  on  the  two-stage  light-gas  guns  at  AEDC  and  the 
specific  problem  to  be  solved. 
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2.  BACKGROUND 


2.1  Two-Stage  Light-Gas  Gun 

Figure  1  shows  the  basic  geometry  and  operation  of  a  two-stage  light-gas  gun.  The  ignition  of  a 
chemical  propellant  drives  a  piston  down  a  pump  tube  containing  a  reservoir  of  a  light  gas  (often 
hydrogen)  contained  by  a  membrane.  The  piston  compresses  die  gas  nearly  isentropically  to  high 
pressures  and  temperatures.  When  die  pressure  exceeds  the  membrane  strength,  the  membrane 
ruptures,  allowing  a  shock  wave  to  accelerate  the  projectile.  As  the  projectile  moves  down  the  barrel, 
the  pressure  applied  at  die  base  of  die  projectile  slowly  decreases,  because  much  of  potential  energy 
initially  present  in  the  pressure  is  converted  into  the  kinetic  energy  of  the  projectile  and  the  gas.  The 
key  to  obtaining  a  very  high  velocity  requires  a  low-molecular-weight  hot  gas,  to  maximize  the  sound 
velocity,  and  the  use  of  long  barrels  to  contain  the  gas  until  most  of  the  energy  is  imparted  to  the 
projectile.  The  survivability  of  the  projectile  depends  on  the  maximum  pressure  applied.  It  is  better  to 
accelerate  with  a  constant,  but  lower  pressure  than  with  a  larger  initial  pressure  that  drops  thereafter. 
To  apply  a  more  constant  pressure  to  the  base  of  the  projectile,  the  piston  may  be  constructed  of  two 
parts:  a  heavy  back  side  and  a  light,  easily  deformable  front  side.  This  front  side,  deforming  as  an 
incompressible  fluid,  tends  to  accelerate  as  it  moves  into  the  constriction,  thus  accelerating  light-gas 
compression  and,  in  turn,  maintaining  a  more  constant  pressure  on  the  base. 

2.2  Geometry 

The  models  and  sabots  come  in  a  variety  of  shapes  and  sizes.  The  "pyro"  model  shown  in  Figure  2  is 
more  or  less  typical.  The  sabot  is  a  lightweight,  nonmetallic  material  (Lexan)  and  the  model  is  usually 
a  graphite-tipped  light  metal,  which  tends  to  ballast  the  center  of  mass  forward  for  better  aerodynamic 
stability.  The  sabot  deforms  plastically  to  provide  a  more  uniform  load  on  the  model,  to  act  as  a  seal 
to  prevent  the  gas  from  leaking  to  the  front  side  of  the  model,  and  to  act  as  a  stabilizer  to  prevent  the 
projectile  from  rotating.  At  the  end  of  the  barrel,  a  barrier  allows  the  model  to  pass  through  but 
prevents  the  sectioned  sabot  from  traveling  further. 

2.3  Axial  Pressure  Loading 

Consider  the  case  of  a  constant  pressure  suddenly  applied  to  the  end  of  a  free-elastic,  homogeneous 
rod.  A  compressive  wave  passes  through  the  rod  at  the  speed  of  sound  in  the  rod.  On  one  side  of  the 
discontinuity  die  stress  is  equal  to  the  applied  pressure,  p,  with  the  material  velocity  equal  to  a  funda¬ 
mental  velocity,  v  =  pl(pc),  where  p  and  c  are  die  rod  density  and  speed  of  sound,  respectively.  This 
wave  reflects  off  die  free  boundary  producing  a  tensile  wave  that  cancels  die  outgoing  compression 
wave.  On  reflection,  the  material  velocity  behind  the  shock  jumps  to  2v.  Each  time  the  wave  reflects 
off  a  boundary,  die  material  velocity  behind  the  shock  ramps  up  one  fundamental  velocity.  This 
ringing  would  occur  indefinitely  in  an  elastic  rod.  In  reality,  material  viscous  damping  and  inhomo¬ 
geneities  will  tend  to  disperse  and  dissipate  the  shock.  Eventually,  the  compressive  stress  field  would 


TR-547/11-91 


3 


Figure  1.  Two-Stage  Accelerated-Reservoir  Light-Gas  Gun 


Figure  2.  Pyro  Model/Sabot  Configuration 
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come  to  a  steady  state,  dropping  linearly  from  the  applied  stress  to  zero  on  the  free  end,  and  the  entire 
rod  would  move  at  the  velocity  of  the  center  of  mass. 

In  this  constant-pressure  elastic  example,  the  rod  is  either  in  compression  or  in  a  zero  state  of  stress. 
Since  most  materials  are  very  strong  in  compression,  one  might  expect  that  die  model/sabot  could 
survive  a  very  large  pressure  pulse.  In  the  gas  gun,  the  applied  pressure  is  far  from  constant,  and  this 
in  itself  can  produce  a  tensile  wave  in  the  rod.  In  Figure  3,  the  typical  pressure  history  on  the  base  of 
the  model  is  predicted  by  the  AEDC  launch  code.  The  pressure  history  exhibits  a  series  of  spikes  that 
are  a  result  of  shock  waves  within  the  gas  reflecting  off  the  back  surface  of  the  sabot  and  off  the 
driving  piston.  Consider  the  elastic  rod  case  but  with  an  applied  pressure  that  decreases  rapidly  with 
time.  As  before,  the  compression  wave  reflects  off  the  free  surface,  returning  as  a  tensile  wave,  but 
now,  because  the  applied  pressure  is  weaker,  the  sum  of  the  incoming  tensile  and  outgoing 
compressive  waves  does  not  cancel  but  may  in  fact  be  in  tension.  These  tensile  waves  are  most 
certainly  generated  when  the  model  undergoes  a  sharp  decompression  as  it  exits  the  barrel  into  a  one- 
atmosphere  environment.  Model  geometry,  more  complicated  than  rods,  can  also  produce  tensile 
stresses  or  stress  concentrations  within  the  body.  A  configuration  constructed  with  different  materials 
will  provide  additional  interfaces  by  which  to  transmit  or  partially  reflect  the  stress  wave. 

Since  the  generation  and  timing  of  shock  waves  within  the  light  gas  depend,  in  part,  on  the  movement 
of  the  model/sabot,  it  is  immediately  clear  that  a  numerical  model  must  account  for  the  coupled 
dynamical  responses  in  both  the  solid  model /sabot  and  the  gas. 

2.4  Strain-Rate-Dependent  Materials 

With  the  model/sabot  subject  to  such  large  stresses,  plastic  deformations  with  large  finite  strains  will 
occur,  especially  in  the  soft  sabot.  Near  the  shock  wave  the  strain  rates  may  be  very  high  (in  the 
elastic  example,  the  region  at  the  shock  interface  has  infinite  strain  rates).  Most  materials,  especially 
plastics  such  as  Lexan,  behave  differently  at  high  strain  rates,  with  stresses  dependent  on  the  strain 
rate.  Accounting  for  this  viscoelastic-viscoplastic  behavior  will  be  important  to  the  characterization  of 
the  state  of  stress  in  the  model. 

From  estimates  of  the  maximum  compressive  stresses  acting  on  the  sabot,  it  appears  that  the  ultimate 
compressive  strength  at  these  high  strain  rates  may  be  2  to  3  times  higher  than  the  static  compressive 
strength.  Instrumented  drop  impact  tests  of  Tardif  and  Marquis  (1963)  indicate  that  both  the  stiffness 
and  the  ultimate  strength  increase  with  strain  rate.  In  these  tests,  a  specimen  sitting  atop  a  stiff  anvil 
bar  is  impact  loaded  by  a  stiff  drop-weight  bar  released  from  a  given  height.  The  nr  an  strain  in  the 
specimen  is  determined  from  integration  of  the  strain  rate,  and  the  strain  rate  is  implied  via  elastic 
wave  relation  involving  the  mean  stress.  This  elastic  relation  introduces  a  degree  of  error,  since  the 
speed  of  sound  is  not  a  constant  but  in  fact  varies  within  the  test.  In  addition,  the  strain  rate  deviates 
considerably  from  the  mean  strain  rate  during  the  course  of  the  test.  Nevertheless,  the  tests  provide  an 
important  measure  of  the  rate  dependence  of  plastics  under  compressive  loads.  Figure  4a  shows  the 
derived  compressive  stress-strain  curve  for  various  mean  strain  rates  (drop  heights)  in  Lexan  (labeled 
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Zelux  in  Tardif  and  Marquis).  The  tests  show  a  dramatic  increase  in  the  dynamic  compressive  strength 
over  the  static  strengths.  The  stiffness  also  appears  to  rise,  but  an  accurate  determination  of  the  slope 
is  difficult  since  the  linear  portion  is  so  small.  Nylon  exhibits  similar  results,  as  seen  in  Figure  4b. 
Lexan  and  nylon  are  in  a  class  of  plastic  materials  in  which  the  3%  total  elongation  yield  strength 
grows  approximately  linear  with  strain  rate. 

2.5  Friction  and  Erosion 

The  wall  friction  will  melt,  bum,  and  erode  the  model/sabot  as  it  moves  down  the  barrel  and  will 
eventually  degrade  the  seal,  allowing  the  gas  to  leak  through  to  the  side  and  front  and  slow  the 
projectile.  The  erosion  of  the  projectile  also  gives  it  more  freedom  to  wobble  or  "ballat"  as  it  travels 
down  the  barrel,  producing  additional  undesirable,  unsymmetrical  side  loads  at  the  contact  points 
between  barrel  and  projectile.  In  fact,  the  inside  diameter  of  the  barrel  at  AEDC  is  tapered  slightly  to 
compensate  for  the  wear  on  the  model/sabot  and  hence  maintains  a  good  seal  and  prevents  ballating 
from  occurring.  The  slope  of  the  taper  has  been  selected  by  measuring  the  amount  of  wear  on  a  tested 
model  in  an  untapered  barrel.  Although  superior  performance  with  the  tapered  barrel  has  been 
observed,  no  other  studies  have  been  conducted  to  determine  if  the  wear  is  over-  or  undercompensated 
in  the  current  taper.  A  numerical  simulation  that  accounts  for  the  friction  and  wear  offers  the  potential 
to  further  refine  the  taper  design. 

During  its  acceleration,  the  compressive  axial  stresses  will  tend  to  expand  the  sabot  into  the  barrel  wall 
to  help  maintain  the  seal  but,  in  turn,  produce  a  higher  friction  stress  and  a  faster  erosion  rate. 
Eventually  the  erosion  will  be  sufficient  that  the  seal  cannot  be  maintained  even  when  the  sabot  is 
under  high  axial  compression.  On  the  other  hand,  if  an  axial  tensile  stress  develops  in  the  sabot,  the 
necking  in  the  radial  direction  can  be  sufficient  to  pull  part  of  the  sabot  away  from  the  barrel  wall. 
One  can  imagine  a  scenario  in  which  the  sabot  is  periodically  constrained  and  unconstrained  as  waves 
of  compressive  stresses  pass  through  the  sabot  and  reflea  off  the  free  surface  as  tensile  stresses.  This 
indicates  that  a  ffiaion  model  must  be  developed  that  recognizes  and  moves  from  three  different 
possible  states:  no  contaa  (free),  sliding  contaa,  and  sticking  contaa. 

2.6  Pressure  Blow-By 

Forces  other  than  axial  pressure  force  also  come  into  play.  The  model/sabot  must  keep  a  seal  to 
prevent  "blow-by,"  the  passage  of  high-pressure  gas  between  the  model/sabot  and  the  barrel  wall.  As 
the  sabot  erodes,  the  space  between  the  sabot  and  the  barrel  wall  becomes  sufficiently  large  to  allow 
the  gas  to  pass  along  to  the  front.  Since  the  pressure  in  front  of  the  sabot  is  subatmospheric,  the  gas 
flow  is  choked  at  the  front  with  a  standing  shock.  Because  of  the  choking,  the  gas  pressure  exerted  on 
the  side  of  the  sabot  will  gradually  decrease  in  moving  to  the  front  of  the  sabot.  It  is  well  known  that 
materials  under  compression  do  not  fail  as  easily  when  under  confining  loads,  and  thus  this  side  gas 
pressure  helps  to  inhibit  the  failure  of  the  sabot,  but  any  possible  turbulence  may  initiate  a  ballating 
motion. 
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2.7  Ballating 

It  is  thought  that  the  sabot  ballats,  or  rocks  back  and  forth,  within  the  bore.  This  ballating  is  felt  to  be 
more  severe  as  the  sabot  erodes  and  is  given  more  freedom  to  rotate.  Although  a  mechanism 
responsible  for  ballating  motion  has  not  been  positively  identified,  several  possible  sources  are: 

(1)  An  unsymmetrical  pressure  loading  on  the  base  of  the  sabot  resulting  from  gas  turbulence. 

(2)  An  unsymmetrical  pressure  loading  on  the  base  of  the  sabot  resulting  from  uneven  diaphragm 
rupture. 

(3)  An  eccentric  center  of  mass. 

(4)  An  unsymmetrical  pressure  loading  on  the  sides  of  the  sabot  resulting  from  leaking  gas 
turbulence. 

(5)  A  slight  curvature  in  the  barrel. 


The  second  is  an  unlikely  source,  since  any  uneven  flow  due  to  rupture  would  have  time  to  relax  to 
symmetry,  especially  in  die  later  stages  of  the  acceleration.  The  third  source  can  be  eliminated  by 
careful  design  of  the  model/sabot.  The  fourth  is  not  likely  because,  as  we  will  show,  the  Reynolds 
numbers  are  no  larger  than  2000  in  the  annulus,  not  large  enough  to  promote  any  strong  turbulent 
pressure  fluctuations.  The  fifth  source  is  also  unlikely  because  the  accelerations  resulting  from  the 
curvature  of  the  barrel  are  small,  even  at  the  maximum  velocities.  Thus,  the  first  seems  to  be  the  most 
likely  source. 

The  Reynolds  number,  based  on  bore  diameter,  peak  temperatures,  and  pressures,  is  about  106,  and 
large,  unsymmetrical  structures  generated  via  any  shock  boundary-layer  interactions  could  be 
sustained.  Each  time  a  shock  reflects  and  moves  away  from  the  base  of  the  sabot  one  might  expect 
large  structures  to  be  generated.  The  magnitude  of  the  pressure  variations  due  to  the  structures  is 
probably  on  the  order  of  the  pressure  drop  of  the  shocks  within  the  bore.  If  the  turbulent  structures 
are  near  the  size  of  the  bore  radius,  then  one  might  expect  significant  moments  could  by  applied  to 
torque  the  sabot. 

Any  more  concrete  analysis  regarding  the  ballating  motion  would  require  a  three-dimensional  model 
along  with  a  Arm  understanding  of  the  statistical  nature  of  die  turbulent  forcing.  This  is  beyond  the 
scope  of  this  Phase  I  effort.  Our  efforts  here  have  concentrated  on  die  analysis  of  axisymmetric  bodies 
undergoing  axisymmetric  loadings. 
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3.  MATHEMATICAL  MODEL 


Each  part  of  the  mathematical  model  will  be  discussed  in  the  subsequent  sections.  Most  of  the 
development  will  be  applicable  to  general  three-dimensional  dynamics,  and  only  when  necessary  will 
the  discussion  be  restricted  to  axisymmetric  dynamics. 

3.1  Kinematics 

In  the  standard  definition,  the  position  of  the  material  point  is  the  sum  of  the  initial  material  position 
plus  the  change  (or  deformation)  from  the  original  position.  Since  the  projectile  is  translating  much 
further  than  the  difference  in  deformation  of  neighboring  points,  this  usual  definition  of  the 
deformation  may  result  in  considerable  round-off  problems,  especially  in  single  precision.  In  an  effort 
to  avoid  these  round-off  errors,  a  slightly  new  definition  for  the  deformation  vector  has  been  adopted. 
In  Figure  5,  the  initial,  current,  and  intermediate  configurations  are  shown.  The  absolute  coordinate 
frame  is  at  the  center  of  mass  of  the  initial  configuration  ft". 

The  deformation  vector  can  be  defined  with  reference  to  an  inertial  frame  that  translates  relative  to  the 
absolute  frame  at  the  velocity  of  the  center  of  mass  of  the  body  at  time  r": 

x(x",  r)  =  x(r")  +  (/  -r")  v(r")  +  x®  +  u"(xa,  r)  m 


Figure  5.  Coordinate  Systems  and  Deformation  Definition 
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The  velocity  of  the  material  point  can  also  be  expressed  as 


The  absolute  acceleration  of  a  material  point  is  then  simply  the  relative  acceleration 


£<x*.  0-^00) 


and  the  Green-Lagrange  strain  is  also  only  an  operation  on  the  deformation  u" 


2e  <„•  = 

*'  '  *;  **  Ac;  a; 


From  the  definition  of  the  center  of  mass 


with  the  fixed  total  mass 


x(0  =  T7  f  jfx(KT,t)dSr 

M  Jq- 

M=  f  pW 
Jo* 


The  center  of  mass  becomes 


x(r)  =  x(r") +  (*-/")  v(r")  + 


OdO0 


and  its  velocity 


v(r)  =  v(r")+-^J  p'^(x',t)dcr 


If  one  defines  a  deformation  u  based  on  a  coordinate  system  that  remains  at  the  center  of  mass 

x(x*.  t)  -  x(r)  +  x®  +  u(x",  t) 

then  die  relationship  between  the  two  different  deformation  systems  is 


u(x",  t)  =  u"(x°,  /)  ~  f  pV(x’,  f)  dQT 
M  Jar 
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with  a  similar  definition  for  the  relative  velocity 


In  the  numerical  scheme  to  be  described  later,  after  each  time  step  the  coordinate  system  is  redefined  to 
be  that  at  die  current  center  of  mass  and  translating  with  die  current  center -of-mass  velocity.  The 
deformation  and  relative  velocity  in  the  new  coordinate  frame  can  then  be  found  from  Equations  (10) 
and  (11). 

3.2  Momentum  Laws 

The  familiar  form  of  die  momentum  equation  is 


d2x  _ 
Pp-.V-O 


3.3  Constitutive  Model 


The  constitutive  model  provides  the  relation  between  the  stress  and  the  strain.  The  discussion  here 
will  focus  on  the  strain-rate-dependent  properties  of  plastics,  such  as  Lexan  and  nylon. 

3.3.1  Elastic-  Viscoplastic 

Following  the  usual  convention,  the  rate  of  change  of  the  stress  tensor  is  related  linearly  to  the  elastic 
part  of  the  strain  rate  tensor 

(i3) 

where  the  strain  rate  is  defined  as  the  sum  of  elastic,  viscoplastic,  and  thermal  components 

**/  +  «/+«/  n4) 


The  thermal  strain  is  defined  as  usual 


£{  =a,r 


The  viscoplastic  strain  will  be  defined  using  the  Bodner-Partom  model  (Bodner,  1984;  Bodner  and 
Partom,  1975): 

M  (16) 

where  r,  is  the  deviatoric  stress  and  A  is  a  function  of  the  second  invariant  of  the  deviatoric  stress  J2\ 


A  =  -St  expf-f^iY— T] 
[  l  2n  \3J2J  J 
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The  hardening  yield  stress  z  is  given  as  a  function  of  the  plastic  work  Wf 

( 

i=Zi  +(Z.-Zi)  exp - L  ,  wf  =  <T,£f 

V  *•  )  (18) 

The  material  constants  zg  and  z,  are  related  to  the  yield  and  the  ultimate  stress,  m  controls  the  rate  of 
work  hardening,  n  controls  the  rate  sensitivity  of  the  material,  and  Da  is  the  maximum  plastic  strain 
rate.  Values  of  these  parameters,  shown  in  Table  1  for  Lex  an,  have  been  chosen  that  reproduce  the 
static  stress  curve  provided  by  GE  Plastics  and  give  a  reasonable  estimate  of  die  ultimate  strength. 
The  corresponding  values  shown  for  aluminum  are  taken  from  Bodnet  (1984)  modified  slightly  for  the 
T6  designation.  The  predicted  stress  as  a  function  of  strain  for  various  strain  rates,  under  uniaxial 
stress  tests,  is  shown  in  Figure  6  for  Lexan.  It  is  evident  from  these  curves  that  the  2%  plastic  yield 


Table  1.  Elastic-Vlscoplaatlc  Material  Parameters 


Material 

n 

■ 

P 

lbm/in3 

GTE 

10*  °F‘ 

ItSWB 

B 

m 

H 

Lexan  500 

D 

m 

m 

17.9 

12.0 

38.0 

1.5 

200 

1000 

Aluminum 

7075-T6 

10.0 

0.32 

0.10 

12.5 

85.5 

115.7 

5 

1000 

Figure  6.  Elastlc-Vlscoplastic  Strain  Rate  Sensitivity 
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stress  varies  approximately  as  a  logarithm  of  the  strain  rate.  Although  the  model  predicts  an  increase 
in  yield,  it  fails  to  reproduce  the  dependence  of  the  stiffness  on  die  strain  rate. 


3.3.2  Viscoelastic 

A  viscoelastic  material  exhibits  a  response  that  is  dependent  on  its  previous  loading  history.  Under 
high  strain  rates,  viscoelastic  materials  show  an  apparent  stiffening.  A  linear  viscoelastic  model 
provides  many  of  the  essential  features  of  this  behavior. 


O  —  EE  +  TfE  (19) 

or 

ct  -E  'e  +  r]  'e  (20) 

Assuming  that  the  strain  rate,  strain,  and  stress  are  initially  zero  and  die  loading  rate  d  is  constant, 
then  the  linear  viscoelastic  law  may  be  solved  to  give  the  strain  and  strain  rate: 


(21) 

(22) 


with  r  =  rf/E  the  viscous  relaxation  time.  As  the  stress  rate  approaches  zero,  the  material  behaves  in 
an  elastic  manner,  but  with  increasing  stress  rate,  the  strain  and  strain  rate  show  an  exponential  decay. 
The  strain  rate  approaches  that  implied  by  the  elastic  limit,  and  the  effective  stiffness 


da 

de 


(23) 


is,  initially,  infinite  and  asymptotically  approaches  the  elastic  modulus.  In  multidimensions,  this  linear 
elastic  law  can  be  expressed  as 

cr4  =  Dykj  +  NvEj  ... 


If  we  further  assume  an  isotropic  material  and  that  Poisson's  ratio  is  unaltered  by  the  viscous  terms, 
then  the  relation  can  be  further  reduced  to 


®T(  =  Djji&j  +  TEj) 

which  is  die  viscoelastic  model  used  in  this  study. 


(25) 
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3.3.3  Viscoelastic- Viscoplastic 

It  is  rather  straightforward  to  combine  the  viscoelastic  and  the  elastic-viscoplastic  strain  laws  into  a 
unified  viscoelastic-viscoplastic  model  stated  simply  as 

+«/-«;)  C6) 

This  model  is  characterized  by  a  total  of  eight  parameters  including  two  elastic  constants,  E  and  v,  the 
viscous  relaxation  time  f,  and  the  five  viscoplastic  parameters  described  earlier.  With  the  best  esti¬ 
mates,  from  the  available  data,  the  material  parameters  for  Lexan  are  shown  in  Table  2.  Note  that  the 
viscoplastic  parameters  are  different  from  those  derived  for  the  elastic-viscoplastic  model  in  Table  1 . 

The  rate  dependence  of  Lexan  is  shown  in  Figure  7  under  constant-stress-rate  uniaxial  loading.  The 
curves  are  labeled  in  terms  of  a  strain  rate  that  is  the  stress  rate  divided  by  the  elastic  modulus,  but  the 
strain  rate  varies  significantly  over  die  duration  of  the  test.  The  prediction  of  the  viscoelastic-visco¬ 
plastic  model  compares  favorably  with  the  results  of  Tardif  and  Marquis  (1963),  as  seen  in  Figure  4a. 
The  Lexan  used  in  the  Tardif  and  Marquis  study  has  a  lower  elastic  modulus  (£  =  200,000  psi)  than 
that  quoted  by  GE  Plastics,  so  a  one-to-one  comparison  is  not  possible.  Figure  8  shows  a  comparison 
of  the  viscoelastic,  elastic-viscoplastic,  and  viscoelastic-viscoplastic  material  modes  at  a  constant  stress 
rate  using  the  identical  set  of  material  parameters  of  Table  2. 

The  stress-strain  curves  of  Tardif  and  Marquis  are  deduced  in  part  with  the  help  of  an  elastic  wave 
model  and  may  suffer  from  this  incorrect  paradigm.  A  more  consistent  approach  to  evaluate  the 
viscoelastic-viscoplastic  material  parameters  would  be  to  compare  and  fit  the  results  of  this  model  to 
various  loading  rate  tests. 

By  accurately  representing  the  geometry  and  boundary  conditions  of  the  material  test  conditions  (such 
as  in  the  drop  tester  or  the  water  cannon  tester),  the  finite  element  model  offers  a  powerful  means  to 
evaluate  the  material  parameters.  The  relaxation  time  may  also  be  related  to  the  Q  damping  factor 
(Vassiliou  et  al.,  1984). 


Table  2.  Viscoelastic-Viscoplastic  Material  Parameters  for  Lexan 


E 

10*  psi 

V 

P 

lbm/in3 

CTE 
10*  °f! 

H 

10s  psi 

h 

103  psi 

n 

m 

Do 

s' 

Q 

0.38 

m 

17.5 

m 

15.0 

25.0 

1.5 

300 

10* 
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3.3.4  Failure 


The  constitutive  model  presented  provides  a  means  to  predict  the  relation  between  stress  and  strain  but 
does  not  define  the  point  at  which  failure  is  expected  to  occur.  Tensile  failure  is  usually  defined  in 
terms  of  ultimate  stresses  or  a  maximum  elongation.  A  more  exacting  measurement  can  be  determined 
from  the  fracture  toughness  and  the  size  of  the  flaws  in  the  material.  Compression  failure  (shear 
failure)  is  generally  defined  by  an  ultimate  compressive  stress.  The  problem  with  these  failure 
definitions  is  that  they  are  strain-rate-dependent  and  generally  only  known  for  static  tests.  It  would  be 
desirable  to  have  a  strength  measurement  that  is  independent  of  the  strain  rate,  but  such  a  measure  is 
not  obviously  clear.  Two  possible  candidates  would  be  plastic  work  or  the  value  of  the  pseudo  yield 
stress  z,  as  defined  in  the  Bodner-Partom  model  in  Equation  (18).  Since  plastic  work  is  a  result  of  the 
integration  of  the  history  of  the  yielding,  it  is  less  likely  to  be  strain-rate-dependent. 

3.4  Boundary  Conditions 

3.4.1  Base  Pressure 

At  the  base  of  the  sabot,  a  pressure  boundary  condition  is  imposed.  The  pressure  is  assumed  to  be  a 
given  function  of  time,  which  must  be  determined  either  experimentally  or  via  another  model,  such  as 
the  AEDC  launch  code. 

3.4.2  Contact  Friction  Model 

The  friction  model  differentiates  between  three  possible  states: 

(1)  No  contact  (free) 

(2)  Sliding  contact 

(3)  Sticking  contact 

In  this  axisymmetric  setting,  the  boundary  conditions  at  x°  =  (r°,  z° )  for  each  state  are  as  given  below. 
No  Contact  (Free) 

The  shear  stress  1  and  normal  stress  <7  acting  on  the  sabot  are  zero,  and  the  deformed  radius  of  the  ma¬ 
terial  point  is  less  than  the  radius  of  the  bore,  R,  which  can  be  a  function  of  z  if  the  barrel  is  tapered. 

r  =  0 ,  <r  =  0 ;  z°,  t)<R  @7) 

Sliding  Contact 

The  deformed  radius  of  the  material  point  is  equal  to  foe  bore  radius,  and  foe  shear  stress  is  equal  in 
sign  to  foe  dynamic  friction  stress  and  points  in  foe  opposite  direction  of  foe  axial  velocity.  The 
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normal  stress  exerted  on  the  sabot  must  be  greater  than  zero. 

Hf,  ie,t)  =  R,  r  =  -sign  (vt)/i/x  ;  a  >  0  (28) 

Sticking  Contact 

The  material  point  is  fixed  to  the  barrel  wall.  The  magnitude  of  die  induced  shear  stress  is  less  than 
the  static  friction,  and  the  normal  stress  is  greater  than  zero. 

r(rA,  z ?  /)  =  !?, z«  at  0  =  z*  ;  ,o>Q  (29) 

Sticking  contact  occurs  only  at  die  initial  stages  of  the  launch  and  quickly  moves  to  sliding  contact 
when  the  shear  stress  induced  by  the  sticking  is  larger  than  the  static  friction  stress. 

3.4.3  Side  Pressure 

3.4.3. 1  Mean  Side-Pressure  Distribution 

As  the  sabot  erodes,  an  annulus  gap  between  die  sabot  and  the  barrel  wall  will  develop,  allowing  the 
gas  to  exert  a  side-pressure  load.  The  approach  taken  here  is  a  greatly  simplified  one,  with  the  prime 
objective  being  to  obtain  an  estimate  of  the  pressure  distribution  along  the  annulus  with  limited 
availability  of  data.  An  idealized  representation  of  the  gap  between  the  projectile  and  the  gun  barrel  is 
sketched  in  Figure  9.  It  is  assumed  that  the  gap  width  is  uniform  circumferentially  around  the 
projectile  and  that  the  gap  has  a  constant  cross-sectional  area  in  the  longitudinal  direction.  The 
Eulerian  coordinate  system  used  has  its  origin  based  at  the  center  of  gravity  of  the  projectile.  It  is 
further  assumed  that  the  projectile  travels  at  a  constant  velocity  V  and  is  not  accelerating. 

Isentropic  Adiabatic  Flow 


v 


Figure  9.  Gas  Flow  In  the  Annulus 


TR-547/11-91 


17 


At  the  base  of  the  projectile,  the  pressure  and  temperature  are  denoted  by  p.  and  Tot  respectively. 
Ahead  of  the  projectile  the  very  low  gun-barrel  pressure  is  denoted  by  ph.  We  analyze  the  gas  motion 
in  the  annulus  using  a  one-dimensional  approach.  The  gas  flow  is  considered  to  be  an  adiabatic  flow 
of  a  perfect  gas,  and  the  only  overall  effect  of  viscosity  is  considered  to  be  a  net  drag  force  on  the  gas. 
The  pressure  differential  across  the  two  ends  of  the  annulus  is  sufficiently  high  that  the  gas  flow  will 
be  choked  at  the  exit  and,  furthermore,  the  flow  within  the  annulus  will  be  completely  subsonic.  The 
pressure  at  the  exit  of  the  annulus  is  determined  by  die  amount  of  frictional  effects  within  the  annulus 
(Shapiro,  1953,  Ch.  6,  Vol.  II)-  The  one-dimensional  analysis  used  here  assumes  that  the  only  effect 
of  the  translating  wall  (see  Figure  9)  is  an  increase  in  the  net  wall  friction. 

In  order  to  obtain  qualitative  estimates  of  the  pressure  variation  within  the  annulus,  two  more 
assumptions  are  necessary.  These  assumptions  are: 

(1)  The  surfaces  of  the  annulus  are  completely  rough. 

(2)  The  friction  factor  is  constant  along  the  annulus. 


The  first  assumption  is  probably  quite  justified  due  to  the  scouring  of  the  gun  barrel  by  the  high-speed 
projectile  and  the  narrow  annulus  width.  It  is  shown  later  that  the  pressure  distribution  along  the 
annulus  is  not  significantly  affected  by  the  choice  of  the  friction  factor,  which  is  an  a  posteriori 
justification  for  the  second  assumption.  With  these  assumptions,  the  analysis  is  identical  to  the 
classical  analysis  of  adiabatic  (constant  enthalpy)  flow  in  a  constant-area  duct.  The  flow  properties 
anywhere  along  the  constant-area  duct  turn  out  to  be  solely  a  function  of  the  locai  Mach  number  of  the 
flow  and  the  conditions  at  the  choking  point  (exit  of  the  duct).  Denoting  the  values  of  velocity  v, 
Mach  number  M,  pressure  p,  stagnation  pressure  pa,  and  density  p  at  the  choking  point  by  a 
superscript  asterisk,  the  local  flow  properties  may  be  expressed  in  a  nondimensional  manner  by  (see 
for  example,  Shapiro,  1953,  Vol.  I). 


With  the  definition 

/J  =  l  +  £— J-M2 

2 


The  required  compressible  relations  are 


(30) 


1  -  M2 

yM 2 


(r+i)M2 

20 


(31) 


p_  =  J_  fr+T 

p  My 


(33) 
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(34) 


t  _  r  +  i 
r "  2p 


p_  j_  nr 

p-~  Mir  +  l 


(35) 


(36) 


T  (37) 

Here,  z  is  the  distance  from  die  local  point  to  the  location  where  choking  occurs.  The  dependence  of 
these  nondimensional  ratios  on  die  local  flow  Mach  number  is  shown  graphically  in  Figure  10.  From 
this  figure  and  an  assumed  friction  factor  /  throughout  the  annulus,  it  is  possible  to  obtain  the  Mach 
number  and  static  pressure  distributions  along  the  annulus.  The  solution  strategy  is  as  follows.  First, 
from  an  assumed  friction  factor  and  die  known  annulus  hydraulic  diameter  (gap)  and  length,  the  factor 
4  flJD  is  determined.  From  this  value  and  Figure  10,  the  Mach  number  at  the  inlet  to  the  annulus  is 
determined.  It  is  assumed  that  die  flow  upstream  of  the  annulus  entrance  point  is  isentropic,  i.e.,  no 
frictional  entrance  effects.  Hence,  knowing  the  entrance  Mach  number  and  the  stagnation  conditions 
upstream  of  the  entrance,  all  flow  quantities  at  the  entrance  to  the  annulus  can  be  determined.  From 
the  now-known  values  of  pjpa'  and  pg  at  the  entrance,  p*  may  be  determined.  The  pressure 
distribution  at  any  point  along  the  annulus  may  then  be  evaluated  by  solving  for  the  local  Mach  number 
using  Equation  (31)  and  inserting  this  value  into  Equation  (32). 

Since  the  flow  is  adiabatic,  the  stagnation  temperature  remains  constant  throughout  the  annulus.  Using 
this  fact  and  the  known  annulus  end  Mach  number  (AT  =  1)  at  the  annulus  exit,  T  may  also  be  deter¬ 
mined  if  Ta  is  provided.  From  p\  AT,  and  T *,  all  flow  quantities  at  the  annulus  exit  may  then  be  de¬ 
termined.  Figure  10  may  now  be  used  to  determine  all  flow  quantities  at  any  point  within  the  annulus. 
For  example,  it  is  known  that  the  Mach  number  within  the  annulus  increases  from  the  inlet  value  to  the 
exit  value. 

Graphs  of  the  Mach  number,  static  pressure,  and  wall  friction  variation  within  the  annulus  are  shown 
in  Figures  11a  through  He  for  several  choices  of  friction  factors.  In  order  to  obtain  these  values,  it 
has  been  assumed  that  the  annulus  is  0.005  in.  wide  and  3  in.  long  and  that  die  upstream  stagnation 
pressure  and  temperature  are  30,000  psi  and  1200  K,  respectively.  These  figures  show  that  die  overall 
pressure  distribution  within  the  annulus  is  not  very  sensitive  to  the  choice  of  friction  factor,  which  is 
an  a  posteriori  justification  for  the  assumption  of  a  constant  friction  factor  throughout  the  annulus. 
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It  is  difficult,  without  resorting  to  numerical  gas  simulations,  to  study  the  question  of  how  the  pressure 
distribution  evolves  in  a  dynamic  sense.  If  the  period  of  the  pressure  variations  is  rather  long  com¬ 
pared  to  the  time  require  for  sound  to  pass  die  length  of  the  sabot  (about  35  /is),  then  the  distribution 
behaves  in  a  quasi-steady  state.  In  this  case,  the  time-dependent  pressure  distribution  follows  from  the 
instantaneous  pressure  at  the  base.  In  situations  in  which  this  is  not  die  case,  a  gross  approximation  is 
used  in  which  the  normalized  distribution  is  the  same  as  that  used  in  the  steady-state  case  but  the  pa  is 
taken  at  a  time  offset  by  the  time  required  for  the  sound  to  reach  the  particular  axial  location. 

3.4.3.2  Annulus  Pressure  Fluctuations 

Based  on  these  compressible  expressions,  the  Reynolds  number  of  the  annulus  flow  ranges  from  700  to 
2000,  and  it  is  anticipated  that  the  flow  is  either  laminar  or  transitional.  Thus,  the  main  cause  of  the 
pressure  fluctuations  within  the  annulus  flow  will  be  due  to  shock  waves  that  reflect  off  the  base  of  the 
projectile.  Knowing  the  strength  of  these  impinging  shock  waves  and  the  appropriate  shock  tube 
formulas,  it  is  possible  to  relate  the  pressure  fluctuation  levels  in  the  annulus  to  those  occurring  at  the 
base  of  the  projectile.  A  more  detailed  treatment  of  the  shock  tube  problem  and  relevant  formulas  may 
be  found  in  Anderson  (1982). 

3.5  Side-Wear  Model 

Due  to  the  very  short  duration  of  the  flight  within  the  barrel,  the  heat  generated  by  friction  at  the 
interface  diffuses  only  slighdy  into  the  interior  of  either  the  sabot  or  the  barrel.  The  diffusion  time 
scales  of  Lexan,  relative  to  the  duration  of  the  test,  are  near  lO7,  and  for  steel  the  ratio  is  much  larger, 
lfr2,  but  still  quite  small.  Since  little  heat  is  dissipated,  much  of  the  friction  energy  goes  into  heating  a 
thin  surface  layer  of  the  sabot  to  its  melting  temperature  and  later  to  promoting  the  burning  of  the 
surface.  One  would  expect  that  after  an  initial  rise  to  the  melting  temperature,  a  pseudo-steady  state 
will  exist  in  which  the  heat  that  is  not  lost  to  dissipation  goes  into  melting  and  ablating  the  sabot.  The 
heat  flows  slightly  into  the  sabot,  and  the  new  cool  portion  of  the  barrel  wall  is  continually  introduced 
to  the  interface.  In  the  coordinate  frame  fixed  to  the  sabot,  this  action  will  produce  a  thin  thermal 
boundary  layer  in  the  barrel  wall.  This  boundary  layer  grows  from  zero  to  some  small  thickness  at  the 
sabot  base.  Assuming  such  a  steady  state,  we  provide  without  proof  the  relation  for  the  rate  of 
ablation  of  the  sabot  rs : 


P,c,(Tm  -To)  +  Pfl,  (38) 

Here  the  subscript  b  indicates  barrel  wall  and  the  s  is  for  sabot.  The  expression  is  derived  from  an 
energy  principal.  The  first  term  in  the  numerator  is  the  heat  energy  flux  generated  by  the  friction,  and 
the  second  term  is  the  energy  flux  dissipated  into  the  barrel.  The  first  term  in  the  denominator  is  the 
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heat  dissipated  into  the  sabot,  and  the  last  term  is  the  heat  removed  by  the  melting  and  ablating.  The 
barrel  velocity  term  derives  from  the  boundary-layer  growth  and  is  given  as 


v*(*)  = 


(39) 


where  v  is  the  speed  of  the  sabot,  x  is  the  distance  from  the  top  of  the  sabot,  and  a*  is  the  thermal 
diffusion  of  the  barrel  wall.  Although  no  data  are  available  for  the  friction  coefficient  at  these  high 
speeds,  Bowden  and  Tabor  (1964)  found  that  at  speeds  of  up  to  1800  fps,  foe  friction  between  nylon 
and  steel  behaved  as 


H  =  1.7  v^4 


A  generalization  of  this  expression  will  be  used  here: 


V0  v£vj 


(40) 


(41) 


This  provides  foe  same  trend  as  that  given  by  Bowden  and  Tabor  but  introduces  a  notion  of  foe  static 
coefficient  of  friction. 

Finally,  in  this  Phase  I  effort  we  have  implemented  only  foe  primary  term  in  foe  rate  of  ablation 
expression: 


(42) 


That  is,  foe  sabot  immediately  begins  to  melt  and  all  foe  heat  generated  by  foe  friction  goes  into 
melting  and  ablating  foe  sabot. 
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4.  FINITE  ELEMENT  MODEL 


4.1  FEAMOD  Overview 

FEAMOD  is  an  inhouse  finite  element  analysis  software  system  developed  over  the  past  ten  years  at 
QUEST  Integrated  (Mueller,  1986).  It  is  intended  as  a  very  flexible  finite  element  analysis  package 
dial  can  be  easily  modified  to  include  new  elements.  More  than  29  element  types  are  available, 
including  elements  for: 

(1)  Isotropic,  orthotropic,  and  anisotropic  elasticity. 

(2)  Elastic/plastic  deformation  with  finite  strain  and  incremental  loading  for  large  deformations  with 
general  yield  functions  and  potentials. 

(3)  Isotropic  and  orthotropic  heat  transfer. 

(4)  Navier-Stokes  laminar  fluid  flow. 

(5)  Porous  media  flow. 

FEAMOD  comes  equipped  with  a  macro  language  to  control  and  tailor  the  solution  algorithm  to 
specific  applications.  Special  macros  are  available  for  time  integration  and  incremental  loading. 
FEAMOD  also  has  translators  for  SUPERTAB,  PATRAN,  and  NISA  for  interactive,  menu-based, 
graphical  pre/post  processing.  The  software  also  includes  state-of-the-art,  preconditioned,  sparse- 
matrix  incore  solvers  suitable  for  both  serial  and  vector  processors  (Oppe  et  al.,  1988).  FEAMOD  has 
been  successfully  applied  to  a  wide  variety  of  engineering  problems  and  will  be  used  as  the  base 
program  in  this  study. 

The  finite  element  method  used  in  FEAMOD  is  that  of  Zienkiewicz  (1977)  and  Bathe  (1982),  and  the 
particular  implementation  found  in  FEAMOD  is  folly  discussed  in  Mueller  (1986).  To  support  the 
further  discussion,  we  describe  here  the  result  of  the  finite  element  discretization  of  the  momentum 
equation: 


r  =  R  -  P  -  CrF  -  Mii  =  0 


(43) 


subject  to  constraints 

p  —  A  —  Cu  -  DF  =  0 

Here  u,  R,  P,  and  F  are,  respectively,  nodal  displacements,  applied  boundary  nodal  forces,  internal 
nodal  forces,  and  nodal  forces  (Lagrange  multipliers)  induced  to  enforce  the  contact  constraints.  C  is 
the  constraint  operator,  and  D  is  a  diagonal  matrix  operating  on  the  nodal  constraint  forces.  A  is  an 
initial  nodal  gap  between  contacting  surfaces,  and  r  and  p  are  the  out-of-balance  nodal  residuals  in  the 
momentum  and  the  constraint  equations.  The  uodai  displacement  u  is  based  on  the  inertial  frame,  as 
discussed  in  Section  3.1. 
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4.2  Numerical  Time  Integration  Method 

FEAMOD  offers  two  choices  for  the  integration  scheme:  the  central  difference  method,  which  is 
explicit  with  a  step  size  stability  limit  of  the  time  required  for  the  sound  wave  to  traverse  the  element, 
or  the  Wilson- 0  method,  which  is  implicit  and  is  absolutely  stable.  An  aluminum  1-in.  slug  with  25 
elements  across  the  slug  will  have  an  explicit  time  step  size  limit  of  near  0.1  /ts  to  remain  stable.  If 
the  slug  is  3  in.  long,  then  die  elastic  wave  will  travel  from  the  base  to  die  top  in  about  15  /is  or  about 
150  steps.  This  is  not  an  inordinately  large  number  of  steps  if  one  is  interested  in  phenomena  in  which 
the  wave  traverses  and  reflects  only  a  few  times.  But,  if  one  is  studying  the  physics  over  the  length  of 
the  shot,  approximately  3  ms,  die  explicit  scheme  will  not  be  efficient.  For  these  reasons,  we  have 
chosen  the  implicit  Wilson- 0  method  for  this  study.  The  drawback  to  the  implicit  scheme  is  that  an 
algebraic  system  of  equations  must  be  inverted.  The  implicit  time  supping  will  require  repeated 
solutions  to  a  linear  system  of  equations,  and  the  cost  of  these  solutions  must  be  minimized.  More 
discussion  of  this  point  will  follow  under  the  linear  solver  section. 

The  Wilson-0  method  assumes  that  the  acceleration  is  a  linear  function  in  the  time  interval  (t,  t  +  OAi). 
With  this  assumption,  the  change  in  the  displacement  and  the  change  in  the  velocity  can  be  evaluated. 
The  scheme  used  here  has  been  modified  to  incorporate  the  constraint  equations  with  the  assumption 
that  the  constraint  forces  also  behave  linearly  over  this  interval.  The  momentum  and  constraint 

equations  are  evaluated  at  time  t  +  6At ,  and  the  change  in  the  displacement  and  the  constraint  forces 

over  this  interval  may  be  solved  via  Newton's  iteration 

K&i*  +  C3F*  =  r(r  +  fifcr)  (45) 

CSu°  +  DdF*  =  p(r  +  0Ar)  ^ 

A 

Here  K  is  the  effective  stiffness  matrix.  The  interpolation  of  the  displacement,  velocity,  and 
acceleration  to  the  new  time  t  +  At  is  described  in  Bathe  (1982),  and  the  nodal  constraint  forces  at  the 
new  time  are  simply  linear  interpolations  of  the  nodal  constraints  at  times  r  and  r  +  flAr. 

4.3  Linear  Solver 

FEAMOD  offers  two  tvpes  of  solvers.  One  is  the  standard  direct  solver  based  on  Gaussian  elimination 
and  on  a  matrix  in  skyline  storage  (Zienkiewicz,  1977).  Since  the  elimination  is  without  pivoting,  it  is 
important  that  any  nodal  constraint  forces  be  arranged  after  the  nodal  displacements  to  help  any  ill 
conditioning  resulting  from  the  zero  entries  in  the  D  matrix  [Equation  (44)].  In  two  dimensions  with 
nxn  quadratic  elements,  it  can  be  shown  that  the  direct  solver  will  require  about  5rt*  operations  to 
eliminate,  4a3  operations  to  forward  and  back  solve,  and  use  about  4n*  words  of  storage.  As  n  grows, 
the  operations  and  storage  rise  dramatically.  In  three  dimensions,  the  problem  is  even  more  severe. 
The  elimination  is  the  most  costly  and  is  required  only  once  for  linear  problems.  For  nonlinear 
problems,  the  computation  of  the  stiffness  matrix  and  its  elimination  may  be  required  as  much  as  once 
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each  time  step.  The  backward  and  forward  elimination  is  carried  out  at  least  once  a  time  step 
regardless  of  linear/nonlinear. 

To  reduce  this  cost,  FEAMOD  offers  an  iterative  solver  based  on  a  preconditioned  conjugate  gradient 
method  (Oppe  et  al.,  1988).  Only  the  nonzeros  of  the  matrix  are  stored  along  with  some  pointing 
arrays.  For  the  same  nxn  quadratic  elements  the  operations  to  solve  a  linear  equation  in  k  iterations 
are  about  SQfci2  and  the  storage  is  about  50/I2  words.  Generally  the  number  of  iterations  is  much  less 
than  h.  This  implies  that,  for  sufficiently  large  problems,  the  iterative  solver  will  be  much  more 
efficient  in  terms  of  both  computer  speed  and  storage.  Experience  shows  this  break-even  point  to  be 
around  n  *  10  or  a  total  of  100  quadratic  elements  for  nonlinear  problems.  There  is,  however,  a 
restriction  on  the  iterative  solver.  The  solver  is  guaranteed  to  converge  only  if  die  matrix  is  symmetric 
positive  definite.  With  the  Lagrange  multipliers  included  in  die  problem,  such  a  condition  is  not  met. 

To  avoid  this  problem,  a  new  conjugate  gradient  solver  that  treats  the  Lagrange  multiplier  terms  has 
been  formulated  and  coded.  The  method  uses  a  preconditioning  similar  to  that  of  the  previous  method 
but  includes  terms  based  on  the  penalty  method  to  handle  the  Lagrange  multiplier.  Since  this  penalty 
is  only  in  the  preconditioner  and  is  not  reflected  in  the  solution  itself,  we  expect  no  stability  problems. 
Indeed,  in  stand-alone  tests,  this  new  solver  appears  to  perform  as  well  on  problems  with  multipliers  as 
the  old  method  does  on  problems  without  the  multipliers.  However,  because  of  time  limitations,  this 
new  solver  has  not  been  implemented  in  the  FEAMOD  structure.  For  this  reason,  the  contact  friction 
element  has  been  used  only  in  conjunction  with  the  direct  solver. 

4.4  Discrete  Constitutive  Model 

The  finite  element  model  must  evaluate  the  stress  at  die  current  time  at  the  Gauss  points.  Applying  a 
second-order-accurate  implicit  scheme  to  the  viscoelastic-viscoplastic  law.  Equation  (26),  the  stress  at 
time  t  +  6At  may  be  written  as 

as/A£'[i«*“ (47) 

This  implicit  algebraic  equation  is  solved  using  a  Newton's  iteration.  The  finite  dement  must  also 
evaluate  the  tangent  stiffness  matrix  at  each  Gauss  point.  Differentiating  Equation  (47)  gives  an 
implicit  linear  equation  for  the  tangent  stiffness: 

dAoi  f  _30  &V  del  dAom 
dA €j  +  6At  j  2  ^  do m  dACj 
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4.5  Boundary  Conditions 

4.5.1  Pnesure  Loading 

The  nodal  boundary  load,  due  to  the  gas  pressure  acting  on  the  sabot  base  or  sides,  is  given  by  the 
usual  inner  product  of  the  pressure  stress  and  the  shape  functions: 


(49) 

The  pressure  acting  on  die  base  is  assumed  homogeneous,  while  on  the  sides  the  pressure  is  taken  from 
Equations  (31)  and  (32).  Because  of  the  spatial  dependence  of  these  side  pressures,  the  inner  products 
are  carried  out  on  an  element  written  explicidy  for  this  effort. 

4.5.2  Friction  Model 

In  die  literature,  die  free/contact  friction  boundary  conditions  are  generally  handled  by  either  a 
Lagrange  multiplier  approach  or  using  a  penalty  method  (stiff  dement  method).  The  Lagrange 
multiplier  approach  offers  a  more  general  treatment  of  the  problem  but  calls  for  the  solution  of 
additional  unknowns.  The  penalty  method  does  not  introduce  any  additional  unknowns,  but  die  penalty 
poorly  conditions  die  matrices  to  be  inverted,  and  a  loss  of  precision  is  expected.  With  single 
precision  arithmetic  and  large  matrices,  die  penalty  method  may  have  potential  stability  problems. 
After  some  initial  experimentation  (see  ELMT12  in  Appendix  A)  with  the  penalty  approach,  we  elected 
to  take  die  more  rigorous  and  general  Lagrange  multiplier  approach  and  accept  the  additional  work  to 
implement  this  approach.  This  decision  also  implied  that  the  current  iterative  solver  would  not  be 
applicable  and  would  require  a  modified  formulation. 

The  free/contact  friction  boundary  conditions  are  enforced  using  a  Lagrange  multiplier  technique  with 
the  constraints  being  die  normal  and  tangential  gaps  between  a  boundary  nodal  point  and  the  barrel 
wall.  The  multipliers  are  die  nodal  forces  required  to  impose  the  boundary  displacement  constraint. 
The  development  follows  Katona  (1983)  except  that  the  constraints  are  dynamic.  With  the  normal  and 
tangential  gaps  defined  for  node  i  as 


s*~  R-r«-uH  (50) 

5*  =  ~u*  (51) 

The  nodal  contributions  to  the  constraint  residuals  and  the  C  and  D  matrices  and  any  explicit  boundary 
loads  P  for  each  contact  state  are  as  given  below. 
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Sticking  contact 

The  Lagrange  multipliers  at  node  i,  FM  and  Fm,  are  both  nonzero  and  defined  by  the  constraints 
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Sliding  contact 

The  Lagrange  multiplier  in  the  normal  direction  is  defined  by  the  normal  constraint.  The  magnitude  of 
the  tangential  constraint  force  is  the  frictional  force  with  a  sign  defined  by  the  local  velocity.  The 
treatment  of  the  friction  is  explicit.  In  other  words,  foe  tangential  forces  of  foe  current  step  are  related 
to  foe  friction  force  based  on  foe  normal  force  of  foe  previous  step. 
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P*  =  P  sign  (v*)/^  -  /^  =  0 

(56) 
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(57) 

Free 

The  normal  and  tangential  forces  are  zero. 

Pm  =  ~F*  =  0 

(58) 
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The  logic  for  moving  from  one  state  to  another  is  shown  in  foe  decision  matrix  of  Table  3.  Some  of 
foe  entries  in  foe  matrix  are  not  self-evident.  The  decision  matrix  does  not  allow  a  state  to  go  directly 
from  a  free  to  a  sticking  contact.  Also,  if  the  body  is  in  sliding  contact  at  the  current  time  but  the 
direction  of  foe  motion  has  changed,  it  is  assumed  that  foe  body  must  go  into  sliding  contact. 
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Table  3.  Free/Contact  Decision  Matrix 


Next  time  step  *+ 1  -* 
Current  time  k  i 

free 

sliding  contact 

sticking  contact 

free 

K>  o 

sliding  contact 

^<0 

F-i0,^t>0 

F**  0,i£uf~lS0 

sticking  contact 

^<0 

As  a  final  note,  the  contact  forces  have  been  normalized  with  a  spring  constant  so  that  they  have  the 
same  units  as  displacement.  In  subsequent  numerical  tests,  a  stability  problem  arose  when  the  sabot 
came  in  recontact  with  the  rigid  barrel  wall.  This  problem  was  effectively  eliminated  by  allowing  the 
barrel  wall  to  have  a  large,  but  not  infinite,  spring  constant. 


4.6  Side-Wear  Model 

The  expression  for  die  rate  of  ablation  involves  die  local  normal  stress  acting  on  the  sabot  due  to  the 
barrel  wall  constraints.  The  equivalent  finite  element  representation  of  Equation  (42)  is 


M/j  =  novty  ds 

with  the  mass  matrix  simply  the  inner  product  of  the  shape  functions 

Mii  =  j^pHtyNj  ds 

Now  the  right-hand  side  can  be  related  to  the  nodal  normal  Lagrange  multiplier  giving 

Mvrj  =  pFMv, 


(61) 


(62) 


(63) 


The  system  is  solved  by  first-order  explicit  time  integration,  and  a  lumped  mass  is  used  to  facilitate  the 
computations.  These  expressions  are  included  in  the  contact  friction  element. 
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5.  EXAMPLE  SIMULATIONS  AND  RESULTS 


To  test  and  demonstrate  the  different  capabilities  of  the  finite  element  model,  two  different  geometries 
and  loading  conditions  are  simulated.  In  the  first  simulation  set,  a  simplified  model  of  die  pyro  model 
and  a  Lex  an  sabot  is  analyzed.  In  the  second  set,  a  Lexan  slug,  undergoing  loadings  similar  to  that  in 
an  impulsive  pressure  test  apparatus  under  development  in  a  parallel  project  (Kolld,  1991).  These  tests 
document  the  evolution  of  the  total  finite  element  model  over  the  course  of  die  study. 

5.1  Aluminum  Pyro  Model/Lexan  Sabot 

The  geometric  model  is  a  simplified  representation  of  die  pyro  model  and  a  Lexan  sabot.  The 
dimensions  and  oudine  of  the  geometric  model  are  the  same,  but  die  pyro  model  is  considered  to  be  a 
single  piece  of  aluminum.  In  reality,  the  pyro  model  is  made  up  of  many  separate  parts  and  materials 
and  contains  some  gaps  between  parts.  We  expect  the  finite  element  model  will  tend  to  underpredict 
the  true  stress  in  the  pyro  model  because  of  the  additional  bulk.  In  all  the  following  cases,  the 
interface  between  the  aluminum  and  the  Lexan  is  frictionless  and  allowed  to  slip.  It  is  also  assumed 
that  the  sabot  remains  in  contact  with  the  pyro  model  and  the  barrel  wall  for  the  duration  of  the  test. 

In  the  initial  simulations,  the  barrel  wall  is  also  assumed  to  be  frictionless,  allowing  the  usual  boundary 
restraint  conditions  offered  with  FEAMOD  to  be  applicable.  The  normal  constraint  between  the  sabot 
and  the  pyro  model  was  handled  by  a  multipoint  constraint  (MPC).  At  the  base  of  this  interface,  the 
MPC  forces  die  axial  displacement  of  nodes  on  either  side  of  the  interface  to  be  the  same.  On  the  side 
interface,  the  MPC  enforces  the  radial  displacement  to  be  the  same.  This  incurs  a  slight  error,  as  the 
side  is  shifted  from  vertical  by  about  10°,  but  simplifies  the  solution  algorithm  because  interfaces 
aligned  with  cylindrical  coordinates  do  not  require  the  solution  of  Lagrange  multiplier  constraints. 

The  model  consisting  of  413  elements  and  1368  nodes  for  a  total  of  2736  degrees  of  freedom  is  shown 
in  Figure  2.  A  fine  grid  is  used  near  die  large  cavity  to  resolve  die  large  stress  gradients.  The  elastic 
model  behaves  properly  if  the  grid  is  coarse  in  these  regions,  but  the  nonlinear  models  fail  to  converge 
if  these  regions  are  not  resolved.  This  will  be  true  of  any  regions  containing  reentrant  corners  or 
cavities.  This  finite  element  model  requires  about  1.2  Mbytes  using  die  iterative  solver  and  5.1 
Mbytes  with  the  direct  solver.  The  generation  of  the  geometry,  nodes,  elements,  and  boundary 
restraints  and  constraints  is  carried  out  with  the  NISA  n  preprocessor.  The  NISA  file  is  then 
translated  into  the  SUPERTAB  universal  file,  which  is  the  standard  FEAMOD  file.  The  FEAMOD 
output  file  is  translated  back  into  a  NISA  file  for  postprocessing.  The  storage  requirements  of  the 
direct  solver  are  dependent  on  the  numbering  scheme.  NISA  does  not  offer  a  utility  for  renumbering 
to  minimize  the  storage,  like  some  packages  such  as  SUPERTAB.  For  this  reason,  the  storage  for  the 
direct  solver  is  higher  than  its  optimal  value. 

The  load  history  on  the  sabot  base  assumes  a  linear  rise  of  pressure  from  0  to  38  ksi  in  about  150  /ts 
and  then  a  slower  linear  drop  to  24  ksi  in  230  /is  for  a  total  simulation  period  of  380  /is.  This  is  the 


TR-547/11-91 


30 


duration  and  magnitude  of  a  typical  spike  produced  by  a  shock  wave  reflecting  off  the  sabot,  as 
predicted  by  the  AEDC  launch  code. 

Initially  the  entire  body  is  at  rest.  The  initial  stresses  within  the  model  are  related  to  the  slight  1-mil 
fit  of  the  sabot  and  die  barrel  wall.  These  initial  stresses  are  predicted  by  first  executing  the  finite  ele¬ 
ment  code  with  a  macro  command  structure  for  static  loading,  representative  of  placing  the  model/ 
sabot  into  die  bore  of  the  barrel. 

5.1.1  Elastic/Frlctlonloss  Barrel 

The  model  predicts  that  both  the  pyro  and  die  sabot  remained  almost  entirely  in  compression  the  entire 
duration  of  the  spike  and  die  maximum  compression  loads  occurred  near  the  peak  loads.  The  von 
Mises  stress  on  an  exaggerated  deformation  plot  is  shown  in  Figure  12.  The  stresses  are  very  high 
near  the  cavity /Lexan/aluminum  junction  and  also  near  the  side  wall/base  of  the  aluminum.  The  Lexan 
is  severely  strained,  especially  near  the  cavity,  with  strains  over  10%.  It  appears  that  for  the  bulk  of 
the  Lexan,  the  von  Mises  stress  is  between  10  and  20  ksi.  This  indicates  that  most  of  the  Lexan  would 
not  yield  if  the  compressive  strength  at  these  strain  rates  (500  s1)  is  near  20  ksi.  As  expected, 
simulations  in  which  the  sabot  is  unconstrained  produce  much  higher  von  Mises  stresses,  but  such  an 
unconstrained  compression  is  not  likely  to  occur.  If  the  sabot  erodes  away,  then  die  gas  pressure 
leaking  into  the  side  will  se.ve  to  constrain  the  model. 

No  significant  tensile  loads  appear  in  the  download.  The  base  of  the  aluminum  model  experiences 
some  tensile  warping,  but  this  is  small,  less  than  10  ksi.  The  period  for  which  an  elastic  wave 
propagates  fiom  the  base  to  the  top  of  the  sabot  is  about  50  /is.  The  time  for  the  elastic  wave  to 
propagate  across  the  aluminum  pyro  model  is  near  15  /is.  These  times  are  rather  short  compared  to  the 
230  /is  drop  in  the  pressure  spike,  and  so  the  model  behaved  almost  quasi-staticaliy. 

The  pressure  spikes  are  a  result  of  shock  waves  impinging  and  reflecting  off  the  base  of  the  sabot. 
Since  these  shock  waves  are  on  the  order  of  a  few  molecules  thick,  there  is  reason  to  believe  that  the 
upsides  of  these  pressure  spikes  are  of  very  short  duration.  Numerical  simulations  of  the  gas  dynamics 
will  inevitably  smear  the  shock.  The  question  is  whether  the  downside,  or  a  portion  of  the  downside, 
is  just  as  steep.  A  careful  examination  of  the  results  of  die  AEDC  launch  code  is  shown  in  Figure  3. 
This  figure  shows  as  little  as  nine  points  in  a  typical  spike,  so  there  may  be  some  numerical  damping 
on  the  downside  (if  each  point  is  a  time  step).  If  the  downsides  are  steeper  than  predicted  by  the  gas 
model,  then  what  effect  will  this  have  on  the  dynamics  of  the  model/sabot?  To  examine  the  relation 
between  tensile  stresses  and  the  ratio  of  die  download  time  to  the  time  for  sound  to  cross  die  model/ 
sabot,  another  elastic  case  was  conducted  in  which  the  duration  of  the  spike  is  cut  in  half  and  by 
fourth.  Significant  tensile  waves  appeared  in  both  simulations  with  the  largest  tensile  stress  in  the 
aluminum.  Figures  13a  and  13b  show  the  axial  stress  contours  at  a  time  in  which  the  tensile  stresses 
peak  in  the  simulation  for  pressure  spike  durations  of  190  and  95  /is.  The  tensile  stress  in  the 
aluminum  is  near  the  maximum  static  tensile  strength.  Viscoelastic  damping  will  reduce  these  loads 


TR-547/11-91 


31 


TR-547/11-91 


32 


TR-547/11-91 


33 


somewhat,  but  the  model  demonstrates  that  the  appearance  of  these  tensile  loads  is  sensitive  to  the 
duration  of  these  spikes. 

5.1.2  Elastic-  Vlscoplastlc/Frictlonless  Barrel 

To  examine  the  effects  of  adding  viscoplasticity  and  geometric  nonlinearities,  the  simulations  are 
repeated  with  die  nonlinear  model.  To  give  an  indication  of  the  dynamic  nature  of  the  stresses,  the 
axial  stresses  at  equal  intervals  of  time  are  shown  on  an  exaggerated  deformation  plot  in  Figures  14a 
through  14n.  Here  we  see  how  in  the  initial  stages  the  stress  field  is  uniform,  but  then  it  becomes 
quite  complicated  as  the  Lexan  and  sabot  interact.  In  Figure  14b  we  see  how  the  stress  wave  has 
propagated  much  foster  in  the  aluminum  than  in  the  Lexan.  The  results  of  the  elastic-viscoplastic  model 
agree  qualitatively  with  the  elastic  runs,  but  the  maximum  stresses  tend  to  be  reduced  as  the  yielding 
redistributes  the  load.  The  contours  of  plastic  work,  shown  in  Figure  15a  through  15k  show  how  the 
failed  region  grows  in  the  vicinity  of  the  cavity.  Most  of  the  yielding  has  occurred  shortly  after  the 
time  of  maximum  pressure.  The  maximum  contour  is  set  to  the  approximate  plastic  work  at  failure  in 
compressive  uniaxial  static  tests  so  that  the  white  region  has  failed.  The  Lexan  is  expected  to  maintain 
its  integrity  because  the  foiled  region  does  not  span  the  diameter. 

5.1.3  Elastlc/Frlctlon  Side  Wall 

In  this  simulation,  the  friction  between  the  barrel  and  the  sabot  and  the  sabot  and  the  pyro  model  is 
considered.  The  coefficient  of  friction  between  the  barrel  and  the  sabot  is  0.5.  This  coefficient  is  a 
severe  test  of  the  robustness  of  the  friction  model,  but  it  is  probably  a  large  overestimate  of  the  actual 
friction.  The  coefficient  of  friction  between  the  sabot  and  the  pyro  model  is  0.2.  The  contact  and 
friction  are  determined  via  one-dimensional  friction  elements  along  the  interfaces,  but  it  is  assumed 
that  the  interfaces  are  always  in  slipping  contact.  Because  no  sticking  contact  is  included,  the  initial 
pressure  applied  at  the  base  is  about  2  ksi  to  overcome  the  friction  of  the  lip  seal.  Unlike  the  previous 
calculations,  the  13-mil  lip  seal  is  considered  in  these  friction  simulations. 

In  Figure  16,  the  von  Mises  stress  contours  at  the  time  of  maximum  pressure  are  shown  in  an 
exaggerated  deformation  plot.  The  friction  tends  to  increase  the  overall  magnitude  of  the  von  Mises 
stresses  in  the  sabot.  The  warping  of  the  sabot  base,  because  of  the  action  of  friction,  is  evident  in  the 
plots.  The  lip  seal  at  the  base  of  the  sabot  appears  to  be  skewed,  but  this  is  an  artifact  of  exaggerating 
the  severe  deformation  of  the  lip.  In  reality,  all  points  at  the  barrel  interface  are  at  the  same  radius. 

This  large  friction  also  affects  dramatically  the  acceleration  of  the  center  of  mass.  The  acceleration  as 
a  function  of  time  is  plotted  in  Figure  17.  The  acceleration  for  a  frictionless  barrel,  which  follows 
from  the  pressure  loading  history,  is  also  shown  for  comparison.  The  friction  acceleration  exhibits 
periodic  oscillations  as  the  compressive  stress  waves  move  through  the  sabot. 
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Figure  14.  Axial  Stress  Contours  (Elastic-Viscoplastic  Model)  (Cont.) 
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Figure  14.  Axial  Stress  Contours  (Elastic- Viscopiastic  Model)  (Cont.) 
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Figure  14.  Axial  Stress  Contours  (Elastlc-Viscoplastic  Model)  (Cont.) 
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Figure  14.  Axial  Stress  Contours  (Elastlc-Viscoplastic  Model)  (Cont.) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastic- Vlscoplastlc  Model) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastlc-Viscoplastlc  Model)  (Cont) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastic-VIscoplaetic  Model)  (Cont.) 


Figure  15.  Plastic  Work  Density  Contours  (Elastlc-Vlscoplastlc  Model)  (Cont.) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastic- Vlscopiastlc  Model)  (Cont.) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastic-Viscopiastic  Model)  (Cont.) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastic-Viscoplastlc  Model)  (Cont.) 
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Figure  15.  Plastic  Work  Density  Contours  (Elastic-VIscoplastlc  Model)  (Cont.) 
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5.1.4  Elastic/Side  Pressure 

This  simulation  is  identical  to  that  described  in  Section  5.1.1  except  that  the  side  loads  are  a  result  of 
the  mean  gas  pressure  in  a  5-mil  gap  between  the  sabot  and  the  barrel  wall  instead  of  the  constraint  of 
the  rigid  wall.  The  von  Mises  stress  at  the  maximum  load  is  shown  in  Figure  18,  which  corresponds 
to  Figure  12  in  the  rigid-wall  constrained  case.  The  side-pressure  loads  have  the  effect  of  decreasing 
the  deviatoric  loads  near  the  sabot  base  but  increasing  the  deviatoric  loads  near  the  top  of  the  sabot  and 
model.  This  can  be  explained  qualitatively  by  examining  the  effect  on  a  slug.  As  one  can  see  from 
Figure  11a,  the  side-pressure  load  axial  drop  is  far  from  linear,  with  little  pressure  drop  occurring  in 
the  first  half  of  the  annulus.  On  the  other  hand,  the  axial  loads  vary  somewhat  linearly,  and  as  a  result 
the  deviatoric  stresses  are  low  at  die  base  and  higher  as  one  moves  to  the  top  of  die  sabot.  In  the  con¬ 
strained  case,  the  side  loads  vary  almost  linearly,  since  they  are  related  to  the  axial  loads  via  Poisson's 
ratio.  Therefore,  the  deviatoric  stress  will  be  highest  where  the  axial  loads  are  largest  at  the  base. 

5.1.5  Elastic/Erosion 

To  validate  the  erosion  model,  a  2.5-in.-diameter,  5-in.-long,  cylindrical  elastic  Lexan  slug  is 
subjected  to  an  idealized  smooth  pressure  history  as  it  travels  down  the  2.498-in.  bore  in  an  untapered 
barrel.  The  duration  of  the  simulation  is  5  ms  with  a  linear  load  peaking  at  40  ksi  at  2  ms  before 
dropping  to  zero.  The  history  of  the  diameter  of  the  slug,  as  a  function  of  axial  coordinate,  is  shown 
in  Figure  19.  To  be  precise,  this  is  die  rebound  diameter  of  an  elastic  Lexan  slug  if  at  that  point  in 
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Figure  19.  Lexan  Slug  Erosion  History 

time  the  pressure  and  the  barrel  wall  are  removed.  The  diameter  curves  are  at  equal  intervals  of  time, 
and  the  location  of  the  slug  down  the  barrel  is  also  labeled.  Several  things  are  evident  from  this  plot: 


Not  surprisingly,  the  model  predicts  that  the  most  wear  occurs  at  the  base  where  die  compressive 
forces  are  largest. 

The  wear  is  not  a  linear  function  of  time  and  reaches  a  limit  diameter  in  which  enough  material 
has  eroded  that  the  axial  compressive  stresses  are  not  sufficient  to  expand  the  slug  into  the  barrel 
wall.  In  the  last  stages,  the  slug  is  primarily  in  ffee  flight.  This  limit  diameter  (0.02  in.  wear  at 
base)  is  probably  somewhat  insensitive  to  die  erosion  rate  since,  if  it  eroded  faster  initially,  it 
must  eventually  slow  down  because  the  compressive  stress  can  expand  out  only  so  much 
material.  The  peak  pressures  and  the  material's  Poisson's  ratio,  on  the  other  hand,  most  likely 
have  a  strong  influence  on  this  limit  diameter. 

The  wear  is  far  from  linear  along  the  barrel  and  most  occurs  in  the  first  8  ft.  A  linear  tapering 
of  the  bore  0.04  in./68  ft  would  not  be  sufficient  to  hold  die  base  in  place  but  would  certainly 
give  rise  to  additional  wear  to  die  front.  Eventually,  die  wear  *.  die  base  of  the  slug  would 
increase  as  it  recontacted  the  barrel. 

The  wear  predicted  by  the  model  does  not  lead  to  the  football  shapes  observed  in  experiments 
(Powell  et  al.,  1985).  The  observed  front  end  wear  is  a  good  clue  that  ballating  does  occur, 
giving  rise  to  the  wear  at  the  extremities. 


TR-547/11-91 


56 


Before  one  can  make  any  firm  conclusions  from  this  wear  model,  more  realistic  geometries  and  loading 
should  be  used  and  the  effects  of  the  taper  and  the  material  parameters  should  be  examined  closely. 

5.2  Water  Cannon  Testa  on  a  Lexan  Slug 

5.2.1  Tsst  Conditions 

A  6-in. -long,  2.5-in. -diameter  Lexan  slug,  which  is  about  1  mil  larger  than  the  bore  of  the  test 
chamber,  is  accelerated  by  impact  with  an  aluminum  piston.  The  mass  and  stiffness  of  the  piston  is 
such  that  the  Lexan  experiences  a  pressure  history  that  has  a  very  sharp  rise  and  then  relaxes  slowly 
until  the  induced  stress  wave  passes  through  the  length  of  the  slug,  reflects,  and  returns  to  the  base. 
At  this  point,  the  Lexan  slug  and  the  piston  separate  and  die  pressure  drops  to  zero.  In  the  test 
apparatus,  the  piston  strikes  a  disk  of  water  held  in  place  by  the  slug  and  a  thin  diaphragm.  The  water 
serves  to  smooth  out  the  impact  on  the  Lexan  and  provides  a  medium  to  measure  die  pressure  with  a 
pressure  transducer.  A  typical  measured  pressure  history  is  shown  in  Figure  20.  This  pressure  history 
is  idealized  as  a  pulse  that  has  a  linear  rise  time  of  20  /is  and  then  drops  linearly  over  a  period  of  130 
/is  for  a  total  duration  of  150  /is.  The  measured  rise  time  is  closer  to  10  /is,  but  the  contact  elements 
behaved  erratically  because  the  grid  is  not  sufficiently  fine  to  resolve  the  short  impulse.  Regardless, 
the  measured  rise  time  may  be  in  error  because  of  die  limited  frequency  response  of  the  transducer. 

The  finite  element  grid  consists  of  250  solid  elements  with  another  25  contact  elements.  There  are  821 
nodes  and  1774  degrees  of  freedom  including  the  constraint  forces.  The  friction  elements  model  the 
full  range  of  possible  states:  free,  sliding  contact,  and  sticking  contact.  Much  of  the  essential 
dynamics  are  represented  here.  In  loading  the  Lexan  slug  into  the  chamber,  a  lubricant  is  used  to 
reduce  the  friction.  Based  on  the  required  load  to  push  the  Lexan  slug  into  the  chamber,  a  coefficient 
of  friction  of  0.1  is  estimated  and  used  in  the  model.  As  before,  the  initial  stresses  in  die  Lexan  are 
determined  by  a  static  simulation. 

5.2.2  Elastic 

Using  the  static  value  of  the  elastic  modulus,  an  elastic  simulation  has  been  conducted  for  over  a  160- 
/is  time  frame  with  a  peak  pressure  of  25  ksi.  In  Figures  21a  through  21h,  a  series  of  snapshots  of  the 
axial  stress  at  different  times  is  shown  on  an  exaggerated  deformation  plot  in  which  the  scale  of  the 
exaggeration  is  held  constant.  The  stress  wave  has  moved  little  by  the  time  of  die  maximum  pressure 
shown  in  Figure  18a.  The  contraction  of  the  slug  is  very  evident  and  some  slight  deviation  from  one 
dimension  is  observed  because  of  the  friction.  At  r  =  80  /is,  shown  in  Figure  21d,  the  wave  reaches 
the  top  of  the  slug  and,  in  the  next  figure,  a  tensile  wave  begins  to  develop  on  reflection.  As  this 
tensile  wave  grows,  the  Lexan  pulls  way  from  die  barrel  wall  with  the  evident  contraction  in  the 
middle.  At  t  *  160  /is  die  Lexan  is  stretched  to  a  maximum  and  the  maximum  tensile  stresses  are  near 
23  ksi  with  tensile  stresses  spanning  the  diameter  as  high  as  13  ksi. 
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Figure  21.  Elastic  Lexan  Slug  Axial  Stress  History  (p^  ■  25  kai) 


TR-547/11-91 


Figure  21.  Elastic  Lexan  Slug  Axial  Stress  History  (p.^  ■  25  ksf)  (Cont.) 


The  garbage  at  the  top  of  the  Lexan  near  die  interface  is  a  result  of  the  slug  returning  to  contact  with 
the  rigid  barrel  wall.  The  local  radial  velocities  at  die  point  of  recontact  are  on  the  order  of  150  fps, 
so  the  impact  is  quite  energetic,  but  since  the  wall  does  not  give,  the  contacting  sabot  bounces  back  off 
in  the  next  time  step.  The  momentum  pushes  the  Lexan  section  back  into  the  wall,  setting  up  the 
oscillations.  The  integration  method  assumes  continuous  accelerations,  but  die  rigid  wall  produces  a 
radical  change  in  the  radial  acceleration.  To  eliminate  this  problem,  the  wall  is  allowed  to  deform 
based  on  a  stiffness  that  is  about  10  times  larger  than  the  sabot.  This  effectively  smooths  the 
accelerations  and  gives  a  more  realistic  model  of  the  barrel  wall. 

Figure  22a  shows  plots  of  the  acceleration  of  the  center  of  mass  along  with  the  acceleration  if  the  slug 
were  frictionless.  The  deviation  from  the  frictionless  case  is  most  pronounced  as  more  of  the  slug  goes 
into  compression  and  expands  against  the  barrel  wall,  producing  more  friction.  In  the  later  stages,  as 
the  Lexan  under  tension  pulls  away  from  the  barrel  wall,  the  acceleration  approaches  the  frictionless 
value.  In  Figure  22b,  the  velocity  of  the  center  of  mass  and  the  point  at  the  base  adjacent  to  the  barrel 
wall  are  plotted.  The  material  point  quickly  rises  to  a  peak  velocity  and  then  falls  back  to  a  value  less 
than  zero  before  rising  again.  If  the  model  is  propagating  the  signal  at  the  correct  speed,  then  the 
wave  reflecting  off  the  base  should  be  occurring  at  the  same  time  that  the  pressure  history  goes  to  zero 
at  150  /is.  The  fact  that  this  elastic  model  is  representing  the  transfer  time  to  within  about  10  /is 
implies  that  the  addition  of  viscoelastic  and  viscoplastic  effects  should  not  greatly  alter  this  time. 

The  axial  stress  at  160  /is  for  a  peak  load  of  40  ksi  is  shown  in  Figure  23  with  the  contour  scales  based 
on  the  range  of  stress  at  this  time.  Here,  the  peaks  are  near  33  ksi  with  the  tensile  stresses  near  21  ksi 
over  the  diameter.  These  tensile  loads  are  quite  high  and,  since  the  Lexan  slug  did  not  fail  at  these 
loads,  either  the  tensile  strengths  are  quite  high  or  this  elastic  simulation,  which  has  no  damping, 
predicts  too  large  of  a  tensile  stress. 

5.2.3  Elastic-Viscoplastic 

Simulations  similar  to  those  in  the  previous  section  have  been  carried  out  using  the  elastic-viscoplastic 
model.  Qualitatively,  the  results  are  similar  to  the  elastic  simulations,  but  the  peaks  stresses  are  not  as 
large  because  of  the  redistribution  to  yielding.  Figure  24  corresponds  to  Figure  23  except  for  the 
change  in  material  model.  The  peak  axial  tensile  stresses  are  about  26  ksi  with  23  ksi  stresses  across 
the  diameter.  The  plastic  work,  shown  in  Figure  25,  indicates  that  most  of  the  failure  yielding  that  has 
occurred  is  associated  with  the  tensile  load.  Although  the  axial  compression  loads  are  higher,  the 
barrel  wall  constrains  the  expansion  to  reduce  the  von  Mises  stresses.  The  tensile  stress  is 
unconstrained,  so  the  deviatoric  stress  is  larger. 

Figures  26  and  27  are  similar  plots  to  Figures  24  and  25  except  that  peak  pressures  are  at  50  ksi. 
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Figure  22.  Acceleration  and  Velocity  of  the  Center  of  Maee  (Elastic  Model) 
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Figure  26.  Elastlc-Vlscoptastlc  Lexan  Slug  Figure  27.  Elaatlc-Viscoplastlc  Lexan  Slug 

Axial  Stress  at  t  ■  160  ps  (p,^  *  50  ksi)  Plastic  Work  at  t  ■  160  pa  (p,^  *  50  ksl) 


5.2.4  Viscoelastic-Viscoplastic 

The  simulation  is  again  repeated  but  with  the  full  viscoelastic-viscoplastic  model.  The  additional 
damping  significantly  effects  the  results  of  these  simulations.  The  time  series  of  axial  stress  contours 
on  die  exaggerated  deformation  plots  in  Figures  28a  through  28h  show  that  the  Lexan's  elastic  rebound 
is  damped,  with  the  result  that  a  tensile  way  appears  only  in  the  last  stages  of  the  deformation  and  is 
considerably  smaller,  peaking  at  around  7  ksi.  This  could  explain  why  the  Lex  an  slug  never  failed, 
even  at  peaks  loads  of  near  45  ksi.  The  speed  at  which  the  stress  waves  move  across  the  slug  also 
seems  to  be  higher,  but  the  reflected  wave  appears  to  have  slowed  somewhat.  The  plot  of  the 
acceleration  and  velocity  of  the  cento:  of  mass  in  Figures  29a  and  29b  appear  almost  identical  to  the 
elastic  case,  indicating  that  the  friction  on  the  side  wall  changed  only  slighdy.  This  is  most  likely  due 
to  the  assumption  the  effective  Poisson's  ratio  of  the  viscoelastic  stiffness  matrix  is  not  changed.  The 
velocity  at  the  base  of  the  sabot  changes  markedly,  showing  the  expected  damped  behavior  and  only  a 
slight  change  in  the  peak-to-peak  period. 

There  is  still  some  question  of  whether  the  timing  of  the  reflection  of  the  tensile  wave  is  consistent 
with  the  duration  of  the  pressure  pulse.  A  definitive  answer  to  this  question  could  be  carried  out  by 
using  FEAMOD  to  model  the  entire  system  of  the  piston,  water,  and  Lexan  slug.  This  is  certainly 
possible  with  the  capabilities  of  FEAMOD,  and  we  hope  to  examine  this  issue  in  the  second  phase  of 
this  effort.  In  addition,  if  the  velocity  or  deformation  history  of  the  top  of  the  sabot  is  measured,  for 
example,  with  a  fiber-optic  interferometer  in  the  brief  time  before  the  slug  goes  sailing,  then  a  good 
measurement  of  the  damping  could  be  obtained  with  the  help  of  the  finite  element  model. 
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Figure  28.  Vlscoelastlc-Viscoplastic  Laxan  Slug  Axial 
Stress  History  (p.^  «  25  ksi) 


Figure  28.  Vlscoelastlc-Vlscoplastic  Lexan  Slug  Axial 
Stress  History  (pM  ■  25  ksl)  (Cont.) 


time  (microseconds) 


a.  Acceleration  (p__„  *25  ksi) 


time  (microseconds) 

b.  Velocity  (p^,  *  25  ksi) 

Figure  29.  Acceleration  and  Velocity  of  the  Center  of  Mass 
(Viscoelastic-Viscoplastic  Model) 
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6.  CONCLUSIONS 


A  finite  dement  model  FEAMOD  has  been  devdoped  and  tested.  This  modd  contains  much  of  the 
essential  physics  for  describing  the  dynamics  of  a  projectile  launched  by  a  two-stage  light-gas  gun.  In 
addition,  FEAMOD  may  be  used  in  conjunction  with  impact  material  tests  to  evaluate  the  material 
parameters  important  to  the  gas  gun  projectile  dynamics.  Capabilities  added  to  FEAMOD  under  this 
effort  include: 

•  A  realistic  constitutive  modd  of  the  viscodastic-viscoplastic  behavior  applicable  to  a  variety  of 
materials,  including  plastics  and  ductile  metals.  The  viscodastic-viscoplastic  modd  accurately 
fits  the  experimental  results  of  high-strain-rate  impact  tests  of  Lexan  and  aluminum. 

•  A  contact  friction  dement  that  differentiates  between  sticking  contact,  sliding  contact,  and  no 
contact  and  comes  with  a  decision  matrix  to  move  from  one  state  to  another.  This  contact 
element  can  be  used  to  describe  contact  forces  between  the  model/sabot  and  the  barrel  wall  or 
between  the  interfaces  within  die  modd  and  die  sabot.  The  contact  friction  element  has  been 
used  to  simulate  a  Lexan  slug  propelled  by  a  water  cannon  tester.  The  model  accurately 
describes  the  radial  expansion  of  die  Lexan  surface  into  the  bore,  which  in  turn  produces  more 
friction.  The  model  also  accounts  for  the  reduction  of  friction  as  the  reflected  tensile  waves 
radially  contract  the  Lexan,  pulling  the  surface  away  from  the  bore. 

•  An  erosion  model  to  account  for  the  ablating  of  the  sabot  surface  at  the  barrel  wall. 

•  A  side-pressure-loading  dement  that  uses  compressible  flow  rdations  to  account  for  the  pressure 
drop  in  any  gas  leaking  around  the  ablated  sides  of  the  model/sabot. 

FEAMOD  comes  with  a  stable  time  integrator  for  the  temporal  dynamics  and  state-of-the-art  linear 
solvers  for  reducing  the  computation  times  and  storage. 

Several  key  issues  identified  in  this  effort  are  as  follows: 

•  The  duration  of  the  spikes  in  the  pressure  history  rdative  to  the  time  required  for  sound  to 
propagate  through  the  solid  projectile  is  an  important  factor  in  determining  whether  tensile 
waves  will  appear  in  the  modd.  In  an  identical  manner,  the  rdative  duration  of  the 
decompression  ramp  as  the  projectile  exits  the  barrd  determines  if  the  projectile  goes  into 
tension.  The  spikes  in  die  pressure  history  are  related  to  gas  shock  waves  impinging  and 
reflecting  off  the  base  of  die  projectile..  For  the  "pyro"  model,  finite  dement  simulations 
indicate  that  tensile  waves  appeared  in  the  aluminum  for  spike  durations  of  less  than  about  0. 1 
ms.  This  duration  is  about  half  the  spike  duration  predicted  by  the  AEDC  gas  dynamics  launch 
code.  Considering  any  possible  numerical  damping  in  these  models,  the  appearance  of  tensile 
waves  in  the  projectile  cannot  be  definhdy  ruled  out.  Careful  numerical  modding  of  the  gas 
dynamics  should  provide  more  definitive  estimates  of  die  spike  and  the  exit  decompression 
durations. 
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•  Regardless  of  the  question  concerning  tensile  waves,  the  finite  element  analysis  predicts 
significant  compression  yielding  at  the  interface  of  the  model  and  the  sabot  near  the  base  of  the 
model.  FEAMOD  quantifies  the  high  stresses  occurring  at  reentrant  corners,  interfaces,  or  gaps 
between  the  differ ent  parts  of  the  model/sabot.  Most  of  the  yielding  is  associated  with  the  peaks 
of  the  pressure  spikes  as  opposed  to  the  average  load  on  die  projectile.  Plastic  work  appears  to 
provide  a  rate-independent  measure  of  the  compression  damage.  Impact  tests,  in  conjunction 
with  the  finite  element  analysis  of  the  tests,  could  provide  estimates  of  the  critical  plastic  work, 
or  the  plastic  work  at  the  ultimate  failure. 

•  The  viscoplasticity  tends  to  smooth  the  peak  stresses  by  redistributing  the  load  near  regions  that 
have  yielded.  The  viscoelasticity  damps  the  stress  wave,  whether  the  material  is  yielding  or  not, 
and  inhibits  the  reflection  of  strong  tensile  waves.  The  degree  of  the  damping  depends  on  the 
viscous  relaxation  time.  This  relaxation  time  is  an  important  parameter  in  describing  the 
behavior  of  Lexan  and  nylon,  and  good  estimates  of  its  magnitude  are  needed. 

•  The  wear  of  the  Lexan  inside  the  barrel  is  neither  a  linear  function  of  position  along  the  barrel 
nor  the  time  after  launch.  This  has  implications  on  the  current  linear  taper  of  the  bore;  more 
studies  are  needed  to  characterize  the  wear  with  more  realistic  pressure  histories  and  actual  bore 
tapers. 

•  Friction  tends  to  increase  the  deviatoric  stresses  and  increase  the  likelihood  of  failure,  but  it  is  a 
secondary  influence  unless  the  coefficient  of  friction  is  above  0.4. 

•  Side-pressure  loads,  as  a  result  of  gas  leaking  around  the  sides  of  the  Lexan,  reduce  the  shear 
stresses  near  the  base  and  tend  to  move  the  maximum  shear  load  closer  to  the  top  of  the  sabot 
and  in  the  model. 

•  The  most  plausible  explanation  of  the  cause  of  projectile  ballating  is  attributed  to  turbulent 
structures  generated  at  the  base  of  the  projectile  as  shock  waves  reflect  from  the  surface.  If 
these  structures  have  length  scales  of  the  size  of  the  bore  radius  and  if  the  pressure  variations 
within  these  structures  are  of  the  same  magnitude  of  the  shock  pressure  rise,  then  large  torques 
can  be  expected.  Furthermore,  if  the  projectile  has  eroded,  then  these  torques  could  produce 
high  yawing  velocities  to  impact  the  bore  wall.  Modeling  the  effects  of  the  barrel  taper  on  the 
erosion  is  important  if  the  designs  to  prevent  possible  ballating  are  to  be  analyzed. 

Several  practical  issues  in  the  use  of  FEAMOD  have  also  been  identified: 

•  Regions  with  rentrant  corners,  gaps,  and  interfaces  must  be  meshed  sufficiently  fine  to  avoid 
numerical  instabilities  associated  with  the  nonlinear  viscoplastic  model. 

•  Since  projectiles  launched  at  the  AEDC  facilities  generally  Lave  a  fairly  complex  geometry,  a 
preprocessor  with  good  "automated  meshing"  capabilities  is  imperative  if  such  analyses  are  to  be 
manpower  effective.  Also,  if  the  analyst  is  to  quickly  view  the  temporal  results,  automation  in 
the  postprocessor  is  important .  It  is  recommended  that  in  Phase  II,  an  evaluation  of  commercial 
pre/postprocessors  be  conducted  to  decide  which  is  most  appropriate  for  AEDC  personnel  and 
the  FEAMOD  software. 
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•  The  iterative  solver  provides  substantial  improvements  over  the  usual  direct  solver  in  terms  of 
memory  and  CPU  time.  The  iterative  solver  currently  in  FEAMOD,  however,  is  not  applicable 
to  the  sliding/sticking  contact  element.  Since  the  advantages  of  the  iterative  solver  will  be 
especially  crucial  for  the  quick  analysis  of  complicated  projectile  geometries,  it  is  recommended 
that  the  a  modified  iterative  solver  be  incorporated  into  FEAMOD  in  Phase  II. 

•  Temporal  elastic  solutions  using  grids  with  near  3000  degrees  of  freedom  over  a  pressure  history 
of  0.4  ms  require  a  little  over  one  megabyte  and  about  30  minutes  on  die  386/387.  Elastic 
simulations  over  the  duration  that  the  projectile  resides  in  the  barrel  would  require  about  8 
hours.  Nonlinear  simulations  would  increase  the  cost  of  a  full  simulation  to  about  24  hours. 
With  complicated  geometries  and  finer  meshes,  the  run  times  will  be  longer,  and  the  analysts 
may  not  be  satisfied  with  this  turnaround.  Therefore,  we  recommend  that  FEAMOD  be  ported 
to  a  computer  with  a  RISC  architecture.  For  example,  the  same  elastic  simulation  over  the  full 
duration  on  an  IBM  RISC  6000  would  require  about  10  minutes!!  Since  FEAMOD  makes  no 
system  calls  and  is  written  entirely  in  FORTRAN,  portability  is  not  a  problem. 

To  conclude,  this  effort  has  met  or  exceeded  all  the  Phase  I  technical  objectives  set  out  in  the  proposal 
and  demonstrates  the  feasibility  of  the  proposed  approach  for  realistic  dynamic  stress  analysis  of 
projectiles  accelerated  by  a  two-stage  light  gas  gun. 
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APPENDIX  A 
FEAMOD  USER'S  GUIDE 


FEAMOD  has  undergone  many  extensions  and  revisions  since  the  last  update  of  the  user's  guide 
(Mueller,  1986),  the  most  notable  being  the  revisions  made  when  it  was  placed  on  the  386/387  and 
with  the  addition  of  the  iterative  solver.  Section  A.l  of  this  appendix  gives  the  1987  FEAMOD  User  's 
Guide.  Section  A.2  gives  the  amendments  to  this  user's  guide. 

A.1  1987  FEAMOD  User's  Guide 

The  FEAMOD  User's  Guide,  as  published  in  Mueller  (1986),  is  provided  on  the  following  pages, 
labeled  pages  78  through  90  of  document  FD-TR-1001.  The  next  section  of  this  appendix  then 
continues  after  page  90. 
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4.  FEAMOD  USER'S  MANUAL 


This  section  serves  as  the  User's  Manual  for  the  FEAMOD  code;  it  is  meant 
to  be  both  a  reference  and  a  guide  for  using  the  code.  Input  formats  and  the 
use  of  the  macro  commands  to  set  up  and  solve  problems  are  explained  here. 

There  are  two  important  aspects  of  running  the  code:  1)  the  size  of  the  M 
array.  The  code  checks  to  make  sure  the  array  is  large  enough  for  the  problem 
by  comparing  the  size  of  the  variable  MAX  in  the  main  program  EXEC  against  the 
markers  used  in  M.  The  code  sets  up  these  array  markers  during  input  as  the 
storage  requirements  for  the  problem  arrays  are  determined  from  input  variables. 
Thus,  it  is  important  that  the  variable  MAX  be  given  the  same  value  as  the 
dimension  of  M.  This  allows  the  code  to  check  the  problem  size  and  to  stop,  if 
-necessary.  2)  The  precision  level  of  variables  in  the  code.  There  is  a 
variable  IPR  set  in  BLOCK  DATA  that  is  the  ratio  of  the  real  variable  size  to 
the  integer  variable  size,  which  is  very  important  for  assigning  addresses  in 
the  M  array.  The  size  of  variables  is  most  easily  determined  by  the  number  of 
bits  or  bytes  in  the  variable  word.  The  variable  IPR  must  be  an  integer  amount. 
Usually,  it  has  a  value  of  2  for  both  single-  and  double-precision  situations; 
however,  the  operator  should  make  sure  the  correct  value  for  IPR  is  being  used. 

The  following  sections  describe  the  format  of  the  input  files  MACRO, 

PROBLEM,  and  MATERIAL. 

4.1  File  Description 

FEAMOD  requires  seven  input /output  files,  of  which  four  are  output  files, 
three  are  files  input  by  the  user,  and  one  is  a  file  set  up  by  the  execution  of 
SUPERTAB.  The  following  sections  describe  these  files. 

4.1.1  MACRO  File  and  Executive  Macro  Commands 

The  MACRO  file  is  the  first  file  read  by  the  code.  The  main  program  begins 
by  reading  the  MACRO  file.  The  first  card  in  the  file  is  the  problem  title 
card  shown  below. 


Problem  Title  Card 

Card  1  FORMAT  (20A4) 

Column  Description  Variable 

1-80  Problem  title  used  in  printout  HEAD  (1-20) 
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The  rest  of  the  MACRO  file  is  used  to  input  the  executive  and  analysis  macros. 
The  analysis  macro  commands  direct  the  solution  steps  to  be  performed  by  the 
code.  They  follow  the  ANALYS  executive  macro  in  the  MACRO  file  and  are 
described  in  Section  4.3  of  this  report. 

The  executive  macros  are  listed  below.  Macro  commands  must  begin  in  column 
one.  Two  F10.0  variables  read  in  are  used  to  specify  which  dump  sets  the  code 
will  read  or  write.  The  INPUT  and  the  ANALYS  executive  macros  have  certain 
automatic  dump  read  or  write  operations.  The  first  F10.0  field  specifies  which 
dump  set  is  to  be  read  by  these  automatic  operations,  the  second  specifies  which 
one  is  to  be  written.  The  INPUT  section  automatically  writes  data  it  has  read 
into  the  dump  set  specified  by  the  executive  macro.  Future  capabilities  of  the 
code  will  allow  the  INPUT  routine  to  read  in  the  specified  dump  set,  modify  the 
data  using  input  macro  commands,  and  write  the  new  set  of  data  to  a  specific 
dump  set. 


Executive  Macros 
(  Format ( 2A4 ,2X ,2F 10.0)  ) 

Macro  Description 

INPUT  Begins  the  data  input  section.  Transfers  program 

control  to  the  INPUT  routine  which  reads  from  the 
PROBLEM  file. 

ANALYS  Begins  the  problem  analysis  section.  Transfers 

program  control  to  PMACR  which  reads  the  analysis 
macros.  These  macros  direct  the  solution  steps  to  be 
performed  by  the  code  and  are  described  in  Section  4.3. 
The  analysis  macros  are  read  from  the  MACRO  file 
immediately  following  the  ANALYS  executive  macro 
command . 

REPORT  Begins  the  report  section.  Transfers  program  control 

to  the  REPORT  routine.  Currently,  this  routine  uses 
only  the  printing  analysis  macros.  These  macros  are 
read  from  the  MACRO  file  immediately  following  the 
REPORT  executive  macro. 

STOP  Terminates  program  operation. 
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4.1.2  PROBLEM  File 


The  PROBLEM  file  is  used  Co  input  data  needed  to  set  up  a  problem.  The 
code  begins  reading  from  this  file  after  either  the  INPUT  executive  macro  or 
the  MESH  analysis  macro  is  executed. 

The  file  begins  with  several  problem  definition  cards  (see  description 
below)  that  are  read  by  the  INPUT  routine.  This  data  is  used  initially  in 
setting  array  storage  sizes  and,  later,  in  running  the  problem.  Input  macro 
commands  follow  the  problem  definition  cards. 


Problem  Definition  Cards 


Card  1  FORMAT  (20A4) 

Column  Description 

1-80  Convenience  title  (may  contain  anything) 

Card  2  FORMAT  (20A4) 

Column  Description 

1-80  Convenience  title  used  to  describe  the 

next  line  of  input;  this  card  is  usually 
"NUMN  PNUME  LNUMA  T  NDM  NDF  NEN  NAD" 


Variable 
TITLE  ( 1-20) 


Variable 
TITLE  (1-20) 


Card  3  FORMAT  (715) 

Column  Description  Variable 

I- 5  Total  number  of  nodes  in  the  problem  NUMNP 

6-10  Number  of  grid  elements  NUMEL 

II- 15  Number  of  element  groups  NUMAT 

16-20  Number  of  spatial  dimensions  (max  *  3)  NDM 

21-25  Number  of  unknowns  per  node  (max  ■  6)  NDF 

26-30  Maximum  number  of  nodes  per  element  NEN 

31-35  Added  size  to  element  matrices  in  excess  of  NAD 

NDF*NEN  (for  nodeless  variables) 

Card  4  FORMAT  (20A4) 

Column  Descr iption  Variable 

1-80  Convenience  title  used  to  describe  the  TITLE  (1-20) 

next  line  of  input;  this  card  is  usually 
"LSYMMLSTATNMASSLMOVE  SECOND" 


Card  5  FORMAT  (2L5,I5,L5) 

Column  Description 

Variable 

1-5 

True,  if  stiffness  matrix  is  symmetric 

LSYMM 

6-10 

True,  if  the  problem  is  static 

LSTAT 

11-15 

Describes  the  mass  matrix 

NMASS 

m 

m 

16-20 

0  for  no  mass  matrix 

1  for  lumped  mass  matrix 

2  for  consistent  mass  matrix 

True,  if  the  X-coordinate  array  is  to 

LM0VE 

21-25 

be  updated  after  each  solution 
(this  is  almost  always  true) 

True  if  the  problem  is  second  order-in  time 

SECOND 
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4.1.3  MATERIAL  File 

The  MATERIAL  file  is  used  to  inpuC  element  material  data  that  is  required 
by  the  MATE  input  macro  command;  its  particular  format  depends  on  the 
particular  element.  See  Section  6  for  the  format  of  each  element. 

4.1.4  UNIVERSAL  File 

The  UNIVERSAL  file  is  an  input  file  that  contains  data  which  defines  the 
grid,  element  types,  loading,  and  boundary  restraints.  This  file  may  be 
generated  via  the  Model  Creation  Module  of  SUPERTAB.  Information  on  the  file 

is  accessed  by  the  input  macro  commands  COOR,  ELEM,  BOUN,  FORC,  and  TEMP. 

4.1.5  DUMP  File 

The  DUMP  file  is  an  unformatted  output  file  that  is  essentially  a  dump  of 
the  contents  in  the  memory  of  a  FEAMOD  run.  A  new  dump  file  is  created  for 
each  analysis  DUMP  macro  command.  The  individual  dump  files  are  distinguished 

by  appending  an  index  to  the  dump  file  name.  The  DUMP  file  is  especially 

useful  for  re-starting  lengthy  calculations. 

4.1.6  REPORT  File 

The  REPORT  file  is  a  formatted  output  file  that  records  the  results  of  the 
analysis  and  any  error  messages.  If  requested,  the  input  data  is  also  recorded 
on  this  file.  The  REPORT  file  is  generated  through  the  analysis  macro  commands 
-  PRTC ,  PRTG ,  PRTK,  PRTM,  PRTR,  PRTV,  AND  REAC. 

4.1.7  ANALYSIS  File 

The  ANALYSIS  file  is  a  formatted  output  file  that  writes  analysis  date  in  a 
universal  file  format  so  that  it  may  be  read  by  SUPERTAB.  The  ANALYSIS  file  is 
generated  through  the  UNIV  analysis  macro  command. 

4.2  Xnpttt  Macro  Commands 

The  input  macro  commands  are  found  in  the  MACRO  file  (see  4.1.1).  The 
subroutine  PMESH  is  called  by  the  INPUT  routine  after  the  problem  definition 
cards  are  read.  Then,  PMESH  reads  the  input  macro  commands,  which  are  listed 
below.  Each  macro  directs  PMESH  to  read  a  data  set. 
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Input  Macro  Commands 
(  Format (A4)  ) 


Macro  Command  Description  of 

(Columns  1-4)  the  Input  Data 


COOR 

ELEM 

BOUN 

FORC 

TEMP 

MATE 

PAGE 
PR  IN 
NOPR 
END 


Grid  coordinate  data 
Grid  element  numbering  data 
Boundary  restraint  data 
Prescribed  nodal  force  data 
Nodal  state  variable  data 

Element  group  data;  directs  program  to  read  input 
information  from  the  IOMATL  file. 

Read  in  carriage  control  variable. 

Print  subsequent  input  data  (default  mode). 

Do  not  print  subsequent  input  data. 

Ends  input  macro  string  (required);  returns  the 
code  to  the  executive  macro  level  in  MACRO. 


Input  macro  commands  can  appear  in  any  order  except  for  END,  which 


must  be 


the  last  command.  COOR,  ELEM,  MATE,  and  END  must  always  be  used,  but  use  of 


the  other  macros  is  optional.  Blank  cards  are  allowed  between  input  macros. 


However,  some  macros  must  be  followed  by  a  specific  set  of  input  cards,  and,  in 
these  cases,  where  a  blank  card  may  have  a  special  meaning.  In  particular, 
COOR,  ELEM,  BOUN,  FORC,  TEMP,  MATE,  and  PAGE  require  input  data  in  specific 


formats.  For  this  reason,  the  data  arrays  for  the  BOUN,  FORC,  and  TEMP  macros 


are  first  zeroed,  in  case  they  are  not  used.  Each  of  these  macros  is  discussed 
in  more  detail  under  its  own  subheading  in  this  section.  The  PRIN,  NOPR,  and 


END  macros  require  no  data. 


4.2.1  COOR  Input  Macro  Command 

The  COOR  input  macro  is  used  to  direct  the  program  to  read  nodal  coordinate 
data  from  the  UNIVERSAL  file.  This  macro  will  cause  the  program  to  search  the 
UNIVERSAL  file  for  a  data  set  *  15  (nodes)  and  read  in  all  data.  For 
axisymmetric  problems,  the  Y  axis  is  the  axis  of  symmetry,  and  X  is  treated  as 
the  radius. 


4.2.2  ELEM  Input  Macro  Command 

The  ELEM  input  macro  is  used  to  direct  the  program  to  read  data  from  the 
UNIVERSAL  file  that  defines  the  grid  element  connectivity,  element  type,  and 
material  type.  This  macro  will  direct  the  program  to  search  the  UNIVERSAL  file 
for  data  set  ■  71  (connectivity). 


FD-TR-1001 


82 


4.2.3  BOUN  Input  Macro  Conanand 

The  BOUN  input  macro  directs  the  program  to  read  the  UNIVERSAL  file  for 
data  defining  the  prescribed  restraints.  It  searches  for  the  data  set  *  89 
(restraints) . 

4.2.4  FORC  Input  Macro  Command 

The  FORC  input  macro  directs  the  program  to  read  data  from  the  UNIVERSAL 
file  that  define  the  nodal  loading  forces.  It  searches  the  UNIVERSAL  file  for 
data  set  ■  90  (loadsets). 

4.2.5  TEMP  Input  Macro  Command 

The  TEMP  input  macro  directs  the  program  to  read  the  UNIVERSAL  file  for 
nodal  state  variable  data  (a  scalar  value).  It  searches  the  file  for  data 
set  *  55  (analysis  data  at  nodes). 

4.2.6  MATE  Input  Macro  Command 

The  MATE  input  macro  is  used  to  input  the  element  routine  number  and  the 
necessary  data  for  each  element  group.  When  the  MATE  macro  command  is  en¬ 
countered,  the  code  begins  reading  from  the  MATERIAL  file. 

The  MATE  macro  needs  a  set  of  input  cards  for  each  element  group.  There 
must  be  no  blank  cards  for  this  macro.  The  first  two  cards  for  each  element 
group  are  read  by  the  PMESH  routine.  The  first  card  is  a  convenience  title 
card  and  the  second  card  contains  the  element  group  number,  MGRP,  and  the 
element  routine  number  for  this  group,  IEL.  Then,  the  element  routine  number, 
IEL,  is  called  to  read  in  the  rest  of  the  input  cards  for  this  element  group. 
The  cards  may  vary  depending  on  the  element  routine  used.  It  is  important  to 
remember  that  the  main  purpose  of  this  call  to  the  element  routine  is  to  fill 
the  C  array. 


MATE  Macro  Data  Cards 


Card  1  Format (20A4) 

Column  Descr iption  Variable 

1-80  Convenience  title  used  to  describe  the  TITL( 1-20) 

next  input  lines.  Usually  this  card  is 
"GROUP  ELnn  OPTIONAL  ELEMENT  TITLE" 
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Card  2  Format (21 5) 

Column  Description  Variable 

1-5  Element  group  number  MGRP 

6-10  Element  routine  number  assigned  to  IEL 

this  element  group 

4.2.7  PAGE  Input  Macro  Command 

The  PAGE  input  macro  is  used  to  read  in  a  vertical  spacing  carriage 
control  variable.  A  default  value  is  provided  (see  below),  so  this  macro  does 
not  have  to  be  used.  If  used,  there  is  only  one  data  card  required,  and  it 
must  immediately  follow  the  macro  command  (blank  cards  may  follow  the  data 
card).  The  card  contains  the  carriage  control  variable,  0,  read  with  the 
format  (Al).  The  variable  is  used  to  control  the  vertical  formatting  of  the 
beginning  of  major  blocks  of  output.  The  values  given  to  0  may  depend  on  the 
compiler  and  printer  used,  but  typical  values  are: 

blank  *  start  on  next  line  (single  spacing) 

0  *  skip  a  line  and  then  start  (double  spacing)  (default) 

1  ■  start  on  first  line  of  next  page 

+  «  start  where  the  printer  pointer  is  currently  located 

(no  advance) 

4.3  Analysis  Macro  Commands 

The  analysis  macro  commands  are  used  to  define  the  solution  steps  for  a 
problem.  They  are  input  on  the  MACRO  file  following  the  ANALYS  executive 
macro  command.  A  complete  list  of  these  macros  appears  below.  The  input 
format  for  these  macros  requires  that  they  begin  in  column  1.  Restrictions  or 
warnings  on  the  use  of  these  variables,  if  any,  for  a  particular  macro  command 
are  described  in  the  list  below.  The  use  of  these  macros  in  the  solution  of 
particular  problems  is  explained  in  the  next  section. 

Analysis  Macro  Commands 
Format  (2A4,  2X,  2F10.0) 

CONR  -  Performs  convergence  check  on  5R. 

CONU  -  Checks  for  1 16  I  I  <  tolerance  ||AU||. 

mg  mg 

DATA  -  Reads  in  a  value  for  THETA  (6) 
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DISP 

DT 

DUMP 

EIGE 

END 

EXCD 

FORM 

FORMA 

IN  IT 

LOOP 

MASS 

MASSV 

MESH 

MOVE 
NEXT 
NOPR 
PR  IN 

PROP  1. 

PRTC 

PRTG 

PRTK 

PRTM 

PRTR 


PRTU 

REAC 

REACA 

SOLV 


Analysis  Macro  Commands  ( Cont . ) 

-  Prints  the  nodal  coordinates,  displacements, 
velocities,  and  accelerations. 

-  Sets  desired  time  step  At. 

•  Dumps  the  storage  for  a  future  restart  or  report. 

-  Solves  for  X  and  U  and  prints  results. 

-  Returns  to  executive  macros. 

-  Solves  for  A. 

-  Forms  -  P. 

-  Forms  R^c  -  P  -  MA. 

-  Checks  Jacobian  of  element  coordinates  for  +  sign. 

-  Begins  loop. 

-  Forms  M. 

-  Forms  M  +  C. 

-  Allows  re-inputting  of  a  portion  of  the  problem 
definition  file. 

-  Updates  the  nodal  coordinates. 

-  Ends  loop. 

-  Sets  flag  to  not  print  macro  commands  (default). 

-  Sets  flag  to  print  the  macro  commands  as  they  occur. 

-  Reads  one  line  of  proportional  load  data  from  below  END 
macro.  (The  number  1  is  input  in  the  first  F10.0  field.) 

-  Prints  control  variables. 

-  Prints  the  variables  stored  at  the  Gauss  points. 

-  Prints  effective  stiffness  matrix. 

-  Prints  effective  mass  matrix. 

-  Prints  contents  of  DR  vector  (states  whether  it  is  the 
out-of-balance  load  vector,  6R,  or  the  displacement 
update  vector,  6). 

-  Prints  the  DU  vector,  AU. 

-  Calculates  and  prints  the  nodal  reactions  for  P 
(use  only  after  CONV  or  CONU  and  before  UPDA) . 

-  ditto  for  P  +  MA 

-  Solves  for  6  and  updates  AU  and  U. 
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Analysis  Macro  Commands  ( Cont . ) 


SSS 

STIF 

STIFA 

STIFCN 

STIFV 

STIFVA 

STRE 


TIME 

TOL 

UPDA 

W1L1 

WIL2 

WIL3 

WIL4 

UNIV  DISP 
UNIV  STRE 
UNIV  REAC 


-  Updates  V  for  steady-state  ice  calculations. 

-  Forms  K. 

-  Forms  K  +  — r-~-  M. 

'  024c2  ’ 

-  «s  ♦  h  i 


-  Forms 


*  *zr 


c. 


-  Forms  K  +  7-7—  C  + 
-  0At  - 


ttu 


-  Calculates  and  prints  stresses  and  strains  at  the  Gauss 
points.  (Note  that  since  calculations  are  done  for  this 
macro,  the  kinematic  variables  and  initial  stresses  must 
be  consistent;  use  before  UPDA  macro). 

-  Increments  time. 

-  Sets  convergence  tolerance. 

-  Updates  the  state  variables  stored  in  the  GP  array. 

-  Prepares  a  consistent  estimate  for  the  kinematic 
variables  at  time  t+0At . 

-  Updates  V  and  A. 

-  Interpolates  kinematic  variables  back  to  time  t+  t. 

-  Makes  new  estimate  for  acceleration  at  time  t+  t-***  t 

-  Writes  displacement  type  variables  to  a  UNIVERSAL  file 

-  Writes  stress  type  variables  to  a  UNIVERSAL  file 

-  Writes  reaction  type  variables  to  a  UNIVERSAL  file 


4.4  Problem-Solving  Methods 

In  this  section,  specific  types  of  problems  are  presented  along  with  the 
analysis  macro  string  required  to  solve  the  problem.  Note  that  these  macro 
strings  are  only  the  necessary  commands  needed  to  solve  the  problem.  Other 
macros,  such  as  parameter-defining  macros,  printout  macros,  and  dump  macros, 
can  be  placed  where  needed  within  these  strings.  Be  sure,  however,  to  heed 
any  restrictions  or  warnings  on  the  use  of  these  macros  given  in  the  complete 
list  of  analysis  macros  in  Section  4.3. 
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4.4.1  Static  Problems 

Static  problems  may  be  solved  using  either  the  Newton  iteration  technique 
or  the  incremental  load  method. 


4.4. 1.1  Newton  Iteration  Technique 

Static  problems  having  equations  of  the  form 


P  -  R. 

-  ~bc 


(4.4.1) 


can  be  solved  using  the  Newton  iteration  technique.  It  is  assumed  that  the 
velocity,  acceleration,  and  mass  tensors  are  unnecessary  for  a  static  problem. 
The  form  of  this  equation  actually  solved  by  Newton's  method  is 

K6  »  Rbc  -  £  ,  (4.4.2) 

where  6  is  the  correction  increment  of  the  displacement  and  P  is  the  nodal 
internal  stress  vector.  The  necessary  macro  string  to  solve  this  type  of 
problem  is 

LOOP  -  Begins  iteration  loop. 

ST IF  -  Forms  K. 

FORM  -  Forms  -  P. 

SOLV  -  Solves  for  5  and  updates  AU  and  U. 

CONU  -  Checks  for  !l£ll  <  tolerance  llAUll. 

NEXT  “  Ends  iteration  loop. 

END  -  Returns  to  executive  macros. 


4.4. 1.2  Incremental  Load  Method 

An  alternate  method  for  solving  static  structural  problems  that  is 
particularly  well  suited  for  plasticity  is  the  incremental  load  method.  The 
program  does  not  iterate  at  each  step  to  achieve  convergence  but  performs  a 
displacement  correction  resulting  from  incrementally  increasing  the  load.  The 
form  of  the  equation  solved  is 

u  -  sbc  -  l  •  <*•*• 

The  necessary  macro  command  string  to  solve  this  problem  by  incrementing  the 
load  is 

PROP  1.  -  Reads  one  line  of  proportional  load  data  which  has 

been  inserted  following  the  END  macro  command. 

(or  using  MESH  macro  after  TIME  macro). 
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LOOP  -  Begins  artificial  time  loop. 

TIME  -  Increments  artificial  time  and  zeros  out  6. 

LOOP  -  Begins  iteration  loop. 

STIF  -  Forms  K. 

FORM  -  Forms  p(R^c  ~  P) ,  where  p  is  the  proportional  load  factor. 

SOLV  -  Solves  for  £  and  updates  AU  and  U. 

CONU  -  Checks  for  ll£ll  <  tolerance  llAUll. 

NEXT  -  Ends  iteration  loop. 

UPDA  -  Updates  state  variables. 

NEXT  -  Ends  artificial  time  loop. 

END  -  Returns  to  executive  macros. 

A  variety  of  loading  histories  can  be  applied  with  the  form  of  the  proportional 
load  factor  expression  presently  available: 

p  *  A^  ♦  A£t  +  A^sin(A^t  +  (4.4.4) 

where  the  A^  coefficients  are  read  from  the  proportional  load  data,  along 
with  the  exponent  L  and  upper  and  lower  limits  on  time. 

4.4.2  Dynamic  Structural  Problems 

Dynamic  structural  problems  can  be  solved  using  either  the  explicit 
central-difference  time  method,  or  the  Wilson-0  implicit  time  integration 
method  coupled  with  Newton's  iteration  technique. 

4.4.2. 1  Explicit  Central-Difference  Time  Integration  Method 

Dynamic  structural  problems  can  be  solved  using  an  explicit  central- 
difference  time  integration  method.  The  form  of  the  equation  solved  in  this 
case  is 

M*  -  Sfcc  -  E  .  .  <‘-4-5> 

where  M  is  the  global  mass  matrix.  Note  that  there  are  limitations  on  the  step 
size  for  this  type  of  problem  (see  Section  2.4.1).  The  necessary  macro  string 
to  solve  this  problem  for  accelerations,  velocities,  and  displacements  is 
MASS  -  Forms  M. 

mm 

LOOP  -  Begins  time-step  loop. 
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FORM  -  Forms  8^  -  P. 

EXCD  *  Solves  for  A.  Updates  all  kinematic  variables. 

NEXT  -  Ends  time-step  loop. 

END  -  Returns  to  executive  macros. 


4. 4. 2. 2  Wilson-9  Method 

Dynamic  structural  problems  having  an  equation  of  the  form 

MA  +  P  -  R.  (4.4.6) 

—  -  -be 


can  be  solved  using  the  Wilson-6  implicit  time  integration  method  coupled 
with  Newton's  iteration  technique.  This  equation  is  solved  for  a  time  step 
that  is  1.4  times  larger  than  the  desired  time  step.  After  getting  the 
solution  at  this  long  time  step,  the  solution  at  the  desired  time  is  found  by 
interpolation  using  Wilson's  assumption  of  linear  variation  of  acceleration 
with  time.  The  form  of  the  equation  that  is  actually  solved  is 


(4.4.7) 


where  «  is  Wilson's  theta  («1.4)  and  6  is  again  the  correction  increment  of 
the  displacement.  The  necessary  macro  string  to  solve  this  type  of  problem  is 
LOOP  -  Begins  time-step  loop. 

WIL1  -  Prepares  a  consistent  estimate  for  the  kinematic 

variables  at  time  t+  At,. 

LOOP  -  Begins  iteration  loop. 

STIFA  -  Forms  K  +  ■-  — ■s  M. 

9  at 

FORMA  -  Forms  R.  -  P  -  MA. 

SOLV  -  Solves  for  *  and  updates  AU  and  U. 

CONU  -  Checks  for  II  ||  tolerance  II  U||. 

WIL2  -  Updates  V  and  A. 

NEXT  -  Ends  iteration  loop. 

WIL3  -  Interpolates  kinematic  variables  back  to  time  t+  t. 

UPDA  -  Updates  the  state  variables  stored  in  the  GP  array. 

WIL4  -  Makes  new  estimate  for  acceleration  at  time  t+At+d  t. 

NEXT  -  Ends  time-step  loop. 

END  -  Returns  to  executive  macros. 
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4.4.3  Evaluation  of  Fundamental  Eigenvalue 

The  determination  of  the  dominant  (smallest)  eigenvalue  of  the  problem 

KU  -  XKU  ,  (4.4.8) 

where  X  is  the  dominant  eigenvalue,  is  determined  by  inverse  iteration  using 
the  following  string  of  macro  commands: 

ST IF  -  Forms  K. 

MASS  -  Forms  M. 

EIGE  -  Solves  for  X  and  U  and  prints  results. 

END  -  Returns  to  executive  macros. 


Following  these  macros,  the  U  vector  contains  the  normalized  eigenvector 
corresponding  to  the  dominant  eigenvalue. 


4.4.4  Dynamic  First-Order  Time  Problems 

Dynamic  first-order  time  problems  having  an  equation  of  the  form 

dU 

B  +  P  *  R, 

~  dt  -  ~DC 


(4.4.9) 


can  be  solved  using  the  Crank-Nicolson  technique.  The  form  of  the  discrete 
equation  solved  is 


-  P 


(4.4.10) 


where  <5  is  the  increment  of  the  displacement  between  time  steps.  The  necessary 

macro  string  to  solve  this  type  of  problem  is 

LOOP  -  Begins  time  step  loop 

STIFCN  -  Forms  -J-B  +  0K 
At- 

FORM  -  Forms  R,  -  P 

-DC 

SOLV  -  Solves  for  <§  and  updates  U 
NEXT  -  Ends  iteration  loop 

UPDA  -  Updates  the  state  variables  stored  in  the  GP  array 

NEXT  -  Ends  time  step  loop 

END  -  Returns  to  executive  macros 

If  the  problem  is  linear,  the  STIFCN  command  can  appear  outside  the  time  loop. 
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A.2  FEAMOD  Amendments 

A.2.1  I/O  Files 


FEAMOD  now  uses.  10  files  that  are  stored  under  a  common  name  (NAME)  with  different  extensions 
in  the  subdirectory  DF\  (Data  File). 

Input  files  (formatted): 


Macro  command  file 

DFNNAME.MAC 

Problem  definition  file 

DFNNAME.PRO 

Grid  definition  file 

DF\NAME.GRD 

Material  definition  file 

DFNNAME.MAT 

Output  files  (formatted): 

Report  output  file 

DF\NAME.RPT 

Analysis  output  file 

DF\NAME.ANL 

Summary  output  file 

DFXNAME.SUM 

Input/output  files  (unformatted) 

Memory  Dump  file 

DFXNAME.DMP 

Nodal  Point  Dump  file 

DF\NAME.NOD 

Gauss  Point  Dump  file 

DF\NAME.GAU 

Most  of  these  files  are  described  in  Mueller  (1986);  the  new  files  are  defined  as  follows: 

Summary  Output  File 

This  file  contains  a  a  brief  summary  of  selected  variables  determined  by  the  code.  Its  output  varies 
depending  on  the  type  of  elements  used.  For  this  project,  the  summary  file  contains  the  history  of 
center  of  mass  variables,  z  position,  velocity,  acceleration,  and  the  position  and  velocity  of  a  selected 
point. 

Nodal  Point  Dump  File 

This  binary  input/output  file  contains  the  nodal  results  for  all  degrees  of  freedom.  It  is  generally  used 
to  transfer  nodal  data  between  different  runs  that  use  the  same  grid  but  apply  different  boundary 
conditions  or  different  problem  types.  In  this  effort,  the  nodal  dump  file  has  been  used  to  transfer 
deformation  and  constraint  forces  from  a  static  run  to  be  used  as  initial  conditions  in  a  dynamic  run. 
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Gauss  Point  Dump  File 


This  binary  input/output  file  contains  Gauss  point  results  for  all  elements  in  the  grid.  It  is  similar  to 
the  nodal  point  dump  file  in  that  it  is  generally  used  to  transfer  Gauss  point  data  between  different  runs 
that  use  the  same  grid  but  apply  different  boundary  conditions  or  problem  types.  In  this  effort,  the 
nodal  dump  file  has  been  used  to  transfer  stress,  strain,  and  other  quantities  from  a  static  run  to  be 
used  as  initial  conditions  in  a  dynamic  run. 

A.2.2  Macro  Commands 

Some  of  the  macro  commands  have  been  modified  and  there  are  several  additional  macro  command 
since  the  documentation  of  Mueller(1986).  The  modified  commands  are: 

SOLV  TOL 

With  the  addition  of  the  iterative  solver,  the  SOLV  command  also  provides  a  relative  tolerance  (TOL) 
that  the  iterative  solver  must  reach  before  the  solution  is  considered  converged. 

CONRABS  TOL 

The  CONR  command  has  a  second  word  extension  of  CONR  that  commands  FEAMOD  to  compute  the 
absolute  norm  of  the  residual.  If  the  absolute  residual  norm  is  less  than  TOL,  then  the  Newton’s 
iteration  is  considered  converged. 

CONUABS  TOL 

The  CONUABS  command  is  identical  to  the  CONRABS  command  except  the  convergence  test  is  on 
the  absolute  normal  of  the  change  in  the  nodal  values  in  the  last  iteration. 

WRTX 

The  WRTX  command  instructs  FEAMOD  to  write  out  data  to  either  the  DFVNAME.NOD  or  the 
DF\NAME.GAU  file.  It  has  several  second  word  extensions: 

WRTXDISP 

Write  all  nodal  data  to  the  DFVNAME.NOD. 

WRTXVEL 

Write  all  nodal  rate  of  change  data  to  the  DFNNAME.NOD. 

WRTXACC 

Write  all  nodal  rate  of  rate  of  change  data  to  the  DF\NAME.NOD. 

WRTXSTRE 

Write  all  Gauss  point  data  to  the  DFNNAME.GAU. 
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REDX 


The  REDX  command  is  identical  to  the  WRTX  command  except  it  instructs  FEAMOD  to  read  data 
from  either  the  DF\NAME.NOD  or  the  DF\NAME.GAU  file.  It  has  identical  second  word 
extensions: 

REDXDISP 

REDXVEL 

REDXACC 

REDXSTRE 

Another  versatile  addition  to  the  FEAMOD  command  library  is  the  WILD  command.  WILD  is  a  wild 
card  command  that  takes  on  its  meaning  depending  on  the  particular  elements  used.  This  command 
instructs  FEAMOD  to  call  a  user-definable  routine  called  WILD.  In  this  effort,  the  wild  card 
command  has  been  used  to  compute  and  update  the  position  and  velocity  of  the  center  of  mass.  Since 
this  requires  integration,  the  WILD  command  instructs  the  elements  to  carry  out  this  integration  with  a 
Gaussian  quadrature. 

PRTG  is  no  longer  a  command  supported  by  FEAMOD.  The  stress  type  variables  (variables  that  are 
stored  at  the  GAUSS  points)  may  be  written  unformatted  with  the  WRTXSTRE  command  or  they  may 
be  interpolated  to  the  nodal  points  and  written  formatted  with  the  UNIVSTRE. 

A.2.3  Element  Library 

The  elements  described  in  the  following  are  elements  developed  specifically  for  this  effort.  Some  are 
modified  versions  of  previous  elements  and  others  are  completely  new  elements. 

ELMT24 

This  is  an  axisymmetric,  linear,  isotropic,  thermo-elastic  infinitesimal  strain  element.  The  element 
accepts  all  two-dimensional  element  types  in  the  FEAMOD  library.  The  data  definition  and  format  of 
the  data  in  the  material  file  is  as  follows: 

card  3  FORMAT(20A4) 

Convenience  title  usually  "YOUNG  POISSON  DENSITY  CTE" 
card  4  (free  format) 

E,  v,p,CTE 
card  5  FORMAT(20A4) 

Convenience  title  usually  "NGAUSS" 
card  6  free  format 

integer  order  of  Gaussian  integration 
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ELMT25 


This  is  a  one-dimensional,  axisymmetric  3-node  element  that  enforces  the  general  contact/free 
constraints  for  contact  with  a  stiff  wall  aligned  in  the  axial  direction.  The  element  also  computes  the 
erosion  rate  based  on  Equation  (42)  to  determine  the  total  radial  nodal  erosion.  The  local  element 
order  is  shown  in  Figure  A.l  with  the  node  number  3  having  the  largest  axial  value.  The  radial 
coordinates  of  the  nodes  need  not  be  equal.  The  data  definition  and  format  of  the  data  in  the  material 
file  is  as  follows: 

card  3  FORMAT(20A4) 

Convenience  tide  usually  "YOUNG  EPS  RADIUS  XMUO  VO  POWER  RHO  H" 
card  4  (free  format) 

E ,  e  ,R  ,  na  ,va  ,  p  ,  H 

where  each  variable  is  defined  as 

£  :  Young's  modulus  of  contacting  body 

£ :  ratio  of  young's  modulus  of  contacting  body  and  rigid  wall 

R:  radius  of  rigid  wall 

Ha  :  static  coefficient  of  friction 

va :  normalizing  velocity  [see  Equation  (41)] 

p  :  density  of  contacting  body 

H  :  Heat  of  ablation 


ELMT28 

This  is  an  axisymmetric  element,  nonlinear  thermo-viscoelastic-viscoplastic  element  for  finite  strains 
and  large  displacements.  The  element  accepts  all  two-dimensional  element  types  in  the  FEAMOD 
library.  The  data  definition  and  format  of  the  data  in  the  material  file  is  as  follows: 

card  3  FORMAT (20A4) 

Convenience  tide  usually  "YOUNG  POISSON  DENSITY  CTE  TRELAX" 
card  4  free  format 
£ ,  v,p,CTE,r 
card  5  FORMAT(20A4) 

Convenience  tide  usually  "ZO  Z1  M  N  D0:’ 

card  6  free  format 

z0  zi  m  n  Da 

where  the  meaning  of  these  viscoplastic  variables  are  described  in  Section  3.3.1 
card  7  FORMAT(20A4) 

Convenience  title  usually  "NGAUSS" 
card  8  free  format 

integer  order  of  Gaussian  integration 
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ELMT12 


This  is  a  one-dimensional  axisymmetric  element,  6-node  element,  that  enforces  sliding  contact  between 
elements  that  are  at  an  arbitrary  orientation.  The  local  node  numbering  of  each  node  is  shown  in 
Figure  A.2.  Opposing  nodes  need  initially  to  be  in  the  general  vicinity  of  each  other.  The  element 
enforces  this  constraint  using  a  penalty  method.  The  penalty  is  defined  by  the  parameter  1/e,  and  £ 
should  be  set  no  smaller  than  10‘3  for  single  precision  calculations.  The  data  definition  and  format  of 
the  data  in  the  material  file  is  as  follows: 

card  3  FORMAT(20A4) 

Convenience  title  usually  "YOUNG  EPS  XMUO" 
card  4  (free  format) 

£,  e,ti 

where  each  variable  is  defined  as 
£  :  Young's  modulus  of  contacting  body 
e :  ratio  of  Young's  modulus  of  contacting  body  and  rigid  wall 
fi:  coefficient  of  friction 


Figure  A.1 


Figure  A.2 
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Universal  Dataset 
Type:  15 
Name:  Nodes 
Status:  Current 

Written  by:  I-DEAS  Level  3  Frame 

Written  by:  I-DEAS  Level  3  Super tab 

Written  by:  1-DEAS  Level  3  Systan 

Read  by:  1-DEAS  Level  3  Super tab 

Read  by:  I-DEAS  Level  3  Systan 

Written  by:  I-DEAS  Level  2.5  Super tab 

Written  by:  1-DEAS  Level  2.5  Systan 

Revision  date:  22-JAN-1986 


RECORD  1: 
FIELD  1: 
FIELD  2: 
FIELD  3: 
FIELD  4: 
FIELDS  5-7: 


FORMAT(4I10, 1P3E13.5) 

NODE  LABEL 

DEFINITION  COORDINATE  SYSTEM  NUMBER 
DISPLACEMENT  COORDINATE  SYSTEM  NUMBER 
COLOR 

3 -DIMENSIONAL  COORDINATES  OF 
NODE  IN  THE  DEFINITION  SYSTEM 


RECORD  1  IS  REPEATED  FOR  EACH  NODE  IN  THE  MODEL. 


EXAMPLE: 


-1 

15 


1 

0.00000E+00 

0 

0 

8 

0 . 00000E+00 

0 . 00000E+00 

2 

5.00000E-02 

0 

0 

8 

5.00000E-01 

0 . 00000E+00 

100 

-4.50000E+00 

-1 


0 


0 


8  1 . 20000E+0 1  1 . 20000E+0 1 


Universal  Dataset 
Type:  55 

Name:  Data  at  nodes 


Status: 

Current 

Written 

by 

I-DEAS 

Level 

3 

Frame 

Written 

by 

I-DEAS 

Level 

3 

Supertab 

Written 

by 

' 1-DEAS 

Level 

3 

Systan 

Written 

by 

1-DEAS 

Level 

3 

Tdas 

Read 

by 

I-DEAS 

Level 

3 

Frame 

Read 

by 

I-DEAS 

Level 

3 

Supertab 

Read 

by 

I-DEAS 

Level 

3 

Systan 

Read 

by 

I-DEAS 

Level 

3 

Tdas 

Written 

by 

I-DEAS 

Level 

2. 

5  Supertab 

Written 

by 

I-DEAS 

Level 

2. 

5  Systan 

Revision  date:  29-JAN-1986 


Systan  and  Tdas  use  a  slightly  different  form  of  this  dataset.  See 
"Systan  Notes"  and  "Tdas  Notes"  below  for  differences. 


RECORD  1: 
FIELD  1: 


FORMAT  (40A2) 
ID  LINE  1 


RECORD  2: 
FIELD  1 


FORMAT  (40A2) 
ID  LINE  2 


RECORD  3: 
FIELD  1: 


FORMAT  (40A2) 
ID  LINE  3 


RECORD  4: 
FIELD  1 


FORMAT  (40A2) 
ID  LINE  4 


RECORD  5: 
FIELD  1 


FORMAT  (40A2) 
ID  LINE  5 


RECORD  6: 


FORMAT  (6110) 


DATA  DEFINITION  PARAMETERS 


FIELD  1 :  MODEL  TYPE 


UNKNOWN 
STRUCTURAL 
HEAT  TRANSFER 
FLUID  FLOW 


FIELD  2:  ANALYSIS  TYPE 


-3: 


4: 


UNKNOWN 
STATIC 
NORMAL  MODE 

COMPLEX  EIGENVALUE  (FIRST  ORDER) 
(COMPLEX  CONJUGATE  PAIRS  NOT 
INCLUDED) 

COMPLEX  EIGENVALUE  (FIRST  ORDER) 
(COMPLEX  CONJUGATE  PAIRS 
INCLUDED) 

TRANSIENT 


5 :  FREQUENCY  RESPONSE 

6:  BUCKLING 

7:  COMPLEX  EIGENVALUE  (SECOND  ORDER) 

FIELD  3:  DATA  CHARACTERISTIC  0:  UNKNOWN 

1 :  SCALAR 

2:  3  DOF  GLOBAL  TRANSLATION 

VECTOR 

3:  6  DOF  GLOBAL  TRANSLATION  & 

ROTATION  VECTOR 
4:  SYMMETRIC  GLOBAL  TENSOR 

5:  GENERAL  GLOBAL  TENSOR 

FIELD  4:  SPECIFIC  DATA  TYPE  0:  UNKNOWN 

1 :  GENERAL 

2:  STRESS 

3:  STRAIN 

4 :  ELEMENT  FORCE 

5 :  TEMPERATURE 

6:  HEAT  FLUX 

7:  STRAIN  ENERGY 

8:  DISPLACEMENT 

9:  REACTION  FORCE 

10:  KINETIC  ENERGY 

1 1 :  VELOCITY 

12:  ACCELERATION 

13:  STRAIN  ENERGY  DENSITY 

14:  KINETIC  ENERGY  DENSITY 

15:  HYDRO-STATIC  PRESSURE 

16:  HEAT  GRADIENT 

17:  CODE  CHECKING  VALUE 

18:  COEFFICIENT  OF  PRESSURE 

FIELD  5:  DATA  TYPE  2:  REAL 

5 :  COMPLEX 

FIELD  6:  NUMBER  OF  DATA  VALUES  PER  NODE  (NDV) 

RECORDS  7  AND  8  ARE  ANALYSIS  TYPE  SPECIFIC 

GENERAL  FORM 

RECORD  7:  FORMAT  (8110) 

FIELD  1 :  NUMBER  OF  INTEGER  DATA  VALUES 

1  <  OR  -  NINT  <  OR  -  10 
FIELD  2:  NUMBER  OF  REAL  DATA  VALUES 

1  <  OR  «  NRVAL  <  OR  •  12 

FIELDS  3-N:  TYPE  SPECIFIC  INTEGER  PARAMETERS 

RECORD  8:  FORMAT  (6E13.5) 

FIELDS  1-N:  TYPE  SPECIFIC  REAL  PARAMETERS 

FOR  ANALYSIS  TYPE  -  0,  UNKNOWN 


RECORD  7 


FIELD  1:1 
FIELD  2:1 
FIELD  3: ID  NUMBER 

RECORD  8: 

FIELD  1:0.0 

FOR  ANALYSIS  TYPE  -  1 ,  STATIC 

RECORD  7: 

FIELD  1:  1 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 

RECORD  8: 

FIELD  1:  0.0 

FOR  ANALYSIS  TYPE  -  2,  NORMAL  MODE 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  4 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  MODE  NUMBER 

RECORD  8: 

FIELD  1:  FREQUENCY  (HERTZ) 

FIELD  2:  MODAL  MASS 

FIELD  3:  MODAL  VISCOUS  DAMPING  RATIO 

FIELD  4:  MODAL  HYSTERETIC  DAMPING  RATIO 

FOR  ANALYSIS  TYPE  -  3,  COMPLEX  EIGENVALUE  (FIRST  ORDER) 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  6 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  MODE  NUMBER 

RECORD  8: 

FIELD  1:  REAL  PART  EIGENVALUE 
FIELD  2:  IMAGINARY  PART  EIGENVALUE 
FIELD  3:  REAL  PART  OF  MODAL  A 
FIELD  4:  IMAGINARY  PART  OF  MODAL  A 
FIELD  5:  REAL  PART  OF  MODAL  B 
FIELD  6:  IMAGINARY  PART  OF  MODAL  B 

FOR  ANALYSIS  TYPE  -  -3,  COMPLEX  EIGENVALUE  (FIRST  ORDER  ) 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  6 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  MODE  NUMBER 


RECORD  8: 

FIELD  1 :  REAL  PART  EIGENVALUE 
FIELD  2:  IMAGINARY  PART  EIGENVALUE 
FIELD  3:  REAL  PART  OF  MODAL  A 
FIELD  4:  IMAGINARY  PART  OF  MODAL  A 
FIELD  5:  REAL  PART  OF  MODAL  B 
FIELD  6:  IMAGINARY  PART  OF  MODAL  B 

FOR  ANALYSIS  TYPE  «  4,  TRANSIENT 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  TIME  STEP  NUMBER 

RECORD  8: 

FIELD  1:  TIME  (SECONDS) 

FOR  ANALYSIS  TYPE  -  5,  FREQUENCY  RESPONSE 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  FREQUENCY  STEP  NUMBER 

RECORD  8: 

FIELD  1:  FREQUENCY  (HERTZ) 

FOR  ANALYSIS  TYPE  -  6,  BUCKLING 

RECORD  7: 

FIELD  1:  1 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 
RECORD  8: 

FIELD  1:  EIGENVALUE 

RECORD  9:  FORMAT  (110) 

FIELD  1 :  NODE  NUMBER 

RECORD  10:  FORMAT  (6E13.5) 

FIELDS  1-N:  DATA  AT  THIS  NODE  (HDV  REAL  OR 

"  COMPLEX  VALUES) 

FOR  ANALYSIS  TYPE  -  7,  COMPLEX  EIGENVALUES  (SECOND  ORDER) 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  6 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  MODE  NUMBER 


RECORD  8: 

FIELD  I:  REAL  PART  EIGENVALUE 
FIELD  2:  IMAGINARY  PART  EIGENVALUE 
FIELD  3:  REAL  PART  OF  MODAL  MASS 
FIELD  4:  IMAGINARY  PART  OF  MODAL  MASS 
FIELD  5:  REAL  PART  OF  MODAL  STIFFNESS 
FIELD  6:  IMAGINARY  PART  OF  MODAL  STIFFNESS 

RECORDS  9  AND  10  ARE  REPEATED  FOR  EACH  NODE. 

NOTES:  1.  ID  LINES  MAY  NOT  BE  BLANK.  IF  NO  INFORMATION 

IS  REQUIRED,  THE  WORD  "NONE"  MUST  APPEAR  IN 
COLUMNS  1-4. 

2.  FOR  COMPLEX  DATA  THERE  WILL  BE  2*NDV  DATA  ITEMS 
AT  EACH  NODE.  THE  ORDER  IS  REAL  PART  FOR  VALUE  1, 
IMAGINARY  PART  FOR  VALUE  1,  ETC. 

3.  THE  ORDER  OF  VALUES  FOR  VARIOUS  DATA 
CHARACTERISTICS  IS: 

3  DOF  GLOBAL  VECTOR:  X,  Y,  Z 


6  DOF  GLOBAL  VECTOR:  X, 

Y,  Z, 

RX, 

RY,  I 

SYMMETRIC  GLOBAL  TENSOR: 

SXX, 

SXY, 

SYY 

sxz, 

SYZ, 

SZZ 

GENERAL  GLOBAL  TENSOR: 

sxx. 

SYX, 

SZX 

SXY, 

SYY, 

SZY 

SXZ, 

SYZ, 

SZZ 

4.  ID  LINE  1  ALWAYS  APPEARS  ON  PLOTS  IN  OUTPUT 
DISPLAY. 

5.  IF  SPECIFIC  DATA  TYPE  IS  "UNKNOWN,"  ID  LINE  2 
IS  DISPLAYED  AS  DATA  TYPE  IN  OUTPUT  DISPLAY. 

6.  TYPICAL  FORTRAN  I/O  STATEMENTS  FOR  THE  DATA 
SECTIONS  ARE: 

READ (LUN,1 000 )NUM 
WRITE 

1000  FORMAT  (1 10) 

READ (LUN, 1010)  (VAL(I) ,I»1 ,NDV) 

WRITE 

1010  FORMAT  (6E13.5) 

WHERE:  NUM  IS  NODE  NUMBER 

VAL  IS  REAL  OR  COMPLEX  DATA 
ARRAY 

NDV  IS  NUMBER  OF  DATA  VALUES 
PER  NODE 

7.  DATA  CHARACTERISTIC  VALUES  IMPLY  THE  FOLLOWING 
VALUES  OF  NDV: 


3:  COMPLEX  EIGENVALUE  (FIRST  ORDER) 

(COMPLEX  CONJUGATE  PAIRS  NOT 
INCLUDED) 

-3:  COMPLEX  EIGENVALUE  (FIRST  ORDER) 

(COMPLEX  CONJUGATE  PAIRS 
INCLUDED) 

7:  COMPLEX  EIGENVALUE  (SECOND  ORDER) 

FIELD  3:  DATA  CHARACTERISTIC  2:  3  DOF  GLOBAL  TRANSLATION 

VECTOR 

3:  6  DOF  GLOBAL  TRANSLATION  & 

ROTATION  VECTOR 

FIELD  4:  SPECIFIC  DATA  TYPE  8:  DISPLACEMENT 

1 1 :  VELOCITY 
12:  ACCELERATION 

Tdas  notes: 

1 .  DATA  ASSOCIATED  WITH  SDRC  MODAL-PLUS  HAS  THE 
FOLLOWING  SPECIAL  FORM  OF  ID  LINE  5. 


FORMAT  (4110) 

FIELD  1:  REFERENCE  COORDINATE  LABEL 
(1-8000) 

FIELD  2:  REFERENCE  COORDINATE  DIRECTION 
1 :  X  DIRECTION 
-1:  -X  DIRECTION 
2:  Y  DIRECTION 
-2:  -Y  DIRECTION 
3:  Z  DIRECTION 
-3:  -Z  DIRECTION 

FIELD  3:  NUMERATOR  SIGNAL  CODE 
1:  DISPLACEMENT 
2:  VELOCITY 
3:  ACCELERATION 
4:  FORCE 

FIELD  4:  DENOMINATOR  SIGNAL  CODE 
1:  DISPLACEMENT 
2:  VELOCITY 
3:  ACCELERATION 
4:  FORCE 

This  plus  records  7  and  8  are  used  to  calculate  Modal  Mass  and 
Residue  with  the  following  specifications. 

ID  line  5 


Read 


-  Modal  Plus  used  the  information  in  ID  line  5  to  scale  the 
residue  amplitude. 

-  Tdas  prompts  the  user  for  the  reference/response  to  be  used 
for  scaling  the  residue  and  ignores  the  reference  in  ID  line 
5. 

Record  8 

Read: 

-  Modal  Plus  zeros  the  value  of  modal  mass  on  the  read. 

-  Tdas  stores  the  value  of  the  modal  mass  directly  in  the  shape 
header. 

Write: 

-  Both  Modal-Plus  and  Tdas  calculate  and  write  the  modal  mass 
based  on  the  parameter  table  (not  the  value) . 

2.  For  Record  6  Field  4-Specific  data  type,  values  0  through  12  are 
as  defined  above.  13  and  IS  through  19  are: 

13:  excitation  force 

15:  pressure 

16:  mass 

\ 7 :  time 

18:  frequency 

1 9 :  rpm 


Universal  Dacaset 
Type:  56 

Name:  Data  on  elements 
Status:  Current 

Written  by:  I-DEAS  Level  3  Frame 

Written  by:  I-DEAS  Level  3  Supertab 

Read  by:  I-DEAS  Level  3  Supertab 

Written  by:  I-DEAS  Level  2.5  Supertab 

Revision  date:  15-JAN-1986 


RECORD  1: 

FORMAT  (40A2) 

FIELD  1: 

ID  LINE  1 

RECORD  2: 

FORMAT  (40A2) 

FIELD  1: 

ID  LINE  2 

RECORD  3: 

FORMAT  (40A2) 

FIELD  1: 

ID  LINE  3 

RECORD  4: 

FORMAT  (40A2) 

FIELD  1: 

ID  LINE  4 

RECORD  5: 

FORMAT  (40A2) 

FIELD  1: 

ID  LINE  5 

RECORD  6: 

FORMAT  (6110) 

DATA  DEFINITION  PARAMETERS 

FIELD  1:  MODEL 

TYPE 

0 

UNKNOWN 

1 

STRUCTURAL 

2 

HEAT  TRANSFER 

3 

FLUID  FLOW 

FIELD  2:  ANALYSIS  TYPE 

0 

UNKNOWN 

1 

STATIC 

2 

NORMAL  MODE 

3 

COMPLEX  EIGENVALUE 

4 

TRANSIENT 

5 

FREQUENCY  RESPONSE 

6 

BUCKLING 

FIELD  3:  DATA 

CHARACTERISTIC 

0 

UNKNOWN 

1 

SCALAR 

2 

3  DOF  GLOBAL  TRANSLATION 
VECTOR 

3 

6  DOF  GLOBAL  TRANSLATION  & 

ROTATION  VECTOR 

4 

SYMMETRIC  GLOBAL  TENSOR 

5 

GENERAL  GLOBAL  TENSOR 

FIELD  4:  SPECIFIC  DATA  TYPE 


0:  UNKNOWN 

1:  UNKNOWN 


FIELD  5:  DATA  TYPE 


2: 

STRESS 

3: 

STRAIN 

4: 

ELEMENT  FORCE 

5: 

TEMPERATURE 

6: 

HEAT  FLUX 

7: 

STRAIN  ENERGY 

8: 

DISPLACEMENT 

9: 

REACTION  FORCE 

10: 

KINETIC  ENERGY 

11: 

VELOCITY 

12: 

ACCELERATION 

13: 

STRAIN  ENERGY  DENSITY 

14: 

KINETIC  ENERGY  DENSITY 

15: 

HYDRO- STATIC  PRESSURE 

16: 

HEAT  GRADIENT 

17: 

CODE  CHECKING  VALUE 

18: 

COEFFICIENT  OF  PRESSURE 

2: 

REAL 

5: 

COMPLEX 

FIELD  6:  NUMBER  OF  DATA  VALUES  FOR  EACH  POSITION  ON  THE 
ELEMENT  (NDV) 

RECORDS  7  AND  8  ARE  ANALYSIS  TYPE  SPECIFIC 
GENERAL  FORM 


RECORD  7: 
FIELD  1: 

FIELD  2: 

FIELDS  3-N: 

RECORD  8: 
FIELDS  1-N: 


FORMAT  (8110) 

NUMBER  OF  INTEGER  DATA  VALUES 
I  <  OR  «  NINT  <  OR  -  10 
NUMBER  OF  REAL  DATA  VALUES 
1  <  OR  -  NRVAL  <  OR  «  12 
TYPE  SPECIFIC  INTEGER  PARAMETERS 

FORMAT  (6E13.5) 

TYPE  SPECIFIC  REAL  PARAMETERS 


FOR  ANALYSIS  TYPE  -  0,  UNKNOWN 


RECORD  7: 
FIELD  1 
FIELD  2 
FIELD  3 


1 

1 

ID  NUMBER 


RECORD  8: 

FIELD  1:  0.0 


FOR  ANALYSIS  TYPE  -  1 ,  STATIC 

RECORD  7: 

FIELD  1:  1 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 


RECORD  8: 

FIELD  4:  0.0 

FOR  ANALYSIS  TYPE  -  2,  NORMAL  MODE 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  3 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  MODE  NUMBER 

RECORD  8: 

FIELD  1:  FREQUENCY  (HERTZ) 

FIELD  2:  MODAL  MASS 
FIELD  3:  MODAL  DAMPING 

FOR  ANALYSIS  TYPE  -  3,  COMPLEX  EIGENVALUE 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  6 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  MODE  NUMBER 

RECORD  8: 

FIELD  1:  REAL  PART  EIGENVALUE 
FIELD  2:  IMAGINARY  PART  EIGENVALUE 
FIELD  3:  REAL  PART  OF  MODAL  A 
FIELD  4:  IMAGINARY  PART  OF  MODAL  A 
FIELD  5:  REAL  PART  OF  MODAL  B 
FIELD  6:  IMAGINARY  PART  OF  MODAL  B 

FOR  ANALYSIS  TYPE  -  4,  TRANSIENT 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  TIME  STEP  NUMBER 

RECORD  8: 

FIELD  1:  TIME  (SECONDS) 

FOR  ANALYSIS  TYPE  «  5,  FREQUENCY  RESPONSE 

RECORD  7: 

FIELD  1:  2 
FIELD  2:  1 

FIELD  3:  LOAD  CASE  NUMBER 
FIELD  4:  FREQUENCY  STEP  NUMBER 

RECORD  8: 

FIELD  1:  FREQUENCY  (HERTZ) 

FOR  ANALYSIS  TYPE  «  6,  BUCKLING 


RECORD  7: 

FIELD  1:  1 
FIELD  2:  I 

FIELD  3:  LOAD  CASE  NUMBER 
RECORD  8: 

FIELD  1:  EIGENVALUE 

RECORD  9:  FORMAT  (2110) 

FIELD  1:  ELEMENT  NUMBER 

FIELD  2:  NUMBER  OF  DATA  VALUES  FOR  THIS 

ELEMENT  (NVAL) 

RECORD  10:  FORMAT(6E13 .5) 

FIELDS  1-N:  DATA  ON  ELEMENT  (NVAL  REAL  OR 

COMPLEX  VALUES) 

RECORDS  9  AND  10  ARE  REPEATED  FOR  ALL  ELEMENTS. 

NOTES:  1.  ID  LINES  MAY  NOT  BE  BLANK.  IF  NO  INFORMATION  IS 

REQUIRED,  THE  WORD  "NONE"  MUST  APPEAR  IN  COLUMNS 
1-4. 

2.  FOR  COMPLEX  DATA  THERE  WILL  BE  2*NVAL  DATA  ITEMS. 
THE  ORDER  IS  REAL  PART  FOR  VALUE  1 ,  IMAGINARY 
PART  FOR  VALUE  1 ,  REAL  PART  FOR  VALUE  2 ,  IMAGINARY 
PART  FOR  VALUE  2,  ETC. 

3.  THE  ORDER  OF  VALUES  FOR  VARIOUS  DATA 
CHARACTERISTICS  IS: 

3  DOF  GLOBAL  VECTOR:  X,  Y,  Z 


6  DOF  GLOBAL  VECTOR:  X, 

Y,  Z, 

RX, 

RY,  RZ 

SYMMETRIC  GLOBAL  TENSOR 

:  SXX, 

SXY, 

SYY, 

SXZ, 

SYZ, 

SZZ 

GENERAL  GLOBAL  TENSOR: 

SXX, 

SYX, 

SZX, 

SXY, 

SYY, 

SZY, 

SXZ, 

SYZ, 

SZZ 

4.  ID  LINE  1  ALWAYS  APPEARS  ON  PLOTS  IN  OUTPUT 
DISPLAY. 

5.  IF  SPECIFIC  DATA  TYPE  IS  "UNKNOWN,"  ID  LINE  2 
IS  DISPLAYED  AS  DATA  TYPE  IN  OUTPUT  DISPLAY. 

6.  TYPICAL  FORTRAN  I/O  STATEMENTS  FOR  THE  DATA 
SECTIONS  ARE: 

READ  (LUN,  1000)  NUM,  NVAL 
WRITE 

1000  FORMAT  (2110) 

READ  (LUN,  1010)  (VAL(I) ,1-1 ,NVAL) 

WRITE 


1010  FORMAT  (6E13.5) 


WHERE:  NUM  IS  ELEMENT  NUMBER 

NVAL  IS  NUMBER  OF  REAL  OR  COMPLEX 
DATA 

VALUES  FOR  THIS  ELEMENT  (MAX  -  90) 
VAL  IS  REAL  OR  COMPLEX  DATA  ARRAY 

7.  DATA  CHARACTERISTIC  VALUES  IMPLY  THE  FOLLOWING 
VALUES  OF  NDV: 

SCALAR:  1 

3  DOF  GLOBAL  VECTOR:  3 
6  DOF  GLOBAL  VECTOR:  6 
SYMMETRIC  GLOBAL  TENSOR:  6 
GENERAL  GLOBAL  TENSOR:  9 

8.  DATA  ON  2D  TYPE  ELEMENTS  MAY  HAVE  MULTIPLE  VALUES 
THROUGH  THE  ELEMENT  THICKNESS.  IN  THESE  CASES, 

NVAL  *  NDV*NPOS  WHERE  NPOS  IS  NUMBER  OF  POSITIONS 
THROUGH  ELEMENT.  NPOS  IS  ALWAYS  I  FOR  SOLIDS. 

THE  ORDER  OF  THE  DATA  IS  NDV  VALUES  FOR  POSITION  1, 
NDV  VALUES  FOR  POSITION  2,  ETC.  THE  ORDER  OF  THE 
NODES  DEFINES  AN  OUTWARD  NORMAL  WHICH  SPECIFIES  THE 
ORDER  FROM  POSITION  1  TO  NPOS. 

9.  ANY  RECORD  WITH  ALL  O.O'S  DATA  ENTRIES  NEED  NOT 
(BUT  MAY)  APPEAR. 

10.  A  DIRECT  RESULT  OF  9  IS  THAT  IF  NO  RECORD  9  &  10 
APPEARS,  ALL  DATA  FOR  THE  DATA  SET  IS  0.0. 

11.  WHEN  NEW  ANALYSIS  TYPES  ARE  ADDED,  RECORD  7  FIELDS 

1  AND  2  ARE  ALWAYS  >  OR  *=  1  WITH  DUMMY  INTERGER  AND 
REAL  ZERO  DATA  IF  DATA  IS  NOT  REQUIRED.  IF  COMPLEX 
DATA  IS  NEEDED,  IT  IS  TREATED  AS  TWO  REAL  NUMBERS, 
REAL  PART  FOLLOWED  BY  IMAGINARY  POINT. 

12.  DATALOADERS  USE  THE  FOLLOWING  ID  LINE  CONTENTION: 

1.  (80A1)  MODEL  IDENTIFICATION 

2.  (80A 1 )  RUN  IDENTIFICATION 

3.  (80A1)  RUN  DATE/TIME 

4.  (80A1 )  LOAD  CASE  NAME 

FOR  STATIC: 

5.  (17H  LOAD  CASE  NUMBER;,  110) 

FOR  NORMAL  MODE: 

5.  (10H  MODE  SAME,  110,  10H  FREQUENCY,  E13.5) 

13.  MAXIMUM  VALUE  FOR  NDV  IS  9. 

MAXIMUM  VALUE  FOR  NVAL  IS  90. 


14.  TYPICAL  FORTRAN  I/O  STATEMENTS  FOR  PROCESSING 
RECORDS  7  AND  8. 

READ  (LUN,  1000)  MINT,  NRVAL,  (1PAR(1) ,1-1 ,N1NT) 
1000  FORMAT  (8110) 

READ  (LUN, 1010)  (NRVAL (I) ,1-1 .NRVAL) 

1010  FORMAT  (6E13.5) 


Universal  Dataset 
Type :  7 1 
Name:  Elements 
Status:  Current 

Written  by:  I-DEAS  Level  3  Frame 

Written  by:  I-DEAS  Level  3  Supertab 

Read  by:  I-DEAS  Level  3  Super tab 

Read  by:  I-DEAS  Level  3  Systan 

Written  by:  I-DEAS  Level  2.5  Supertab 

Revision  date:  15-JAN-1986 


RECORD  1:  FORMAT (71 10) 

FIELD  1:  ELEMENT  LABEL 

FIELD  2:  FE  GRAPHICAL  DESCRIPTION  ID 

FIELD  3:  FE  DESCRIPTOR  ID 

FIELD  4:  PHYSICAL  PROPERTY  TABLE  NUMBER 

FIELD  5:  MATERIAL  PROPERTY  TABLE  NUMBER 

FIELD  6:  COLOR 

FIELD  7:  NUMBER  OF  NODES  ON  ELEMENT 

RECORD  2:  FORMAT  (8110) 

FIELDS  I-N:  NODE  LABELS  DEFINING  ELEMENT 

RECORDS  1  AND  2  ARE  REPEATED  FOR  EACH  ELEMENT  IN  THE  MODEL. 

EXAMPLE: 

-1 

71 


1 

19 

1 

1 

1 

8 

8 

11 

12 

13 

16 

21 

20 

19 

15 

2 

12 

2 

2 

1 

8 

16 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

• 

46 

• 

• 

124 

19 

1 

1 

1 

8 

8 

9 

10 

11 

15 

19 

18 

17 

14 

-1 


Universal  Dataset 
Type:  88 

Name:  Constraint  set 
Status:  Current 

Written  by:  I-DEAS  Level  3  Supertab 
Read  by:  I-DEAS  Level  3  Supertab 
Revision  date:  15-JAN-1986 


RECORD  1:  FORMAT  (2110) 

FIELD  I:  CONSTRAINT  SET  NUMBER 

FIELD  2:  CONSTRAINT  TYPE 

1  -  KINEMATICS 

2  -  COUPLED  DOFS 

3  -  MULTI-POINT  CONSTRAINT  EQUATIONS 

RECORD  2:  FORMAT  (20A2) 

FIELD  1:  CONSTRAINT  SET  NAME 

FOR  CONSTRAINT  TYPE  -  1  -  KINEMATICS 

RECORD  3:  FORMAT  (2110,612) 

FIELD  1:  KINEMATIC  NODE  LABEL 

FIELD  2:  COLOR  NUMBER 

FIELDS  3-8:  SWITCHES  FOR  DOFS  1-6 

0  -  OFF 

1  -  ON 

RECORD  3  IS  REPEATED  FOR  EACH  KINEMATIC  OF  THE  CONSTRAINT  SET. 

FOR  CONSTRAINT  TYPE  -  2  -  COUPLED  DOFS 

RECORD  3:  FORMAT  (3110,612) 

FIELD  1:  INDEPENDENT  NODE  LABEL 

FIELD  2:  COLOR  NUMBER 

FIELD  3:  NUMBER  OF  DEPENDENT  NODES 

FIELDS  4-9:  SWITCHES  FOR  "Y)FS  1-6 

0  -  OFF 

1  -  ON 

RECORD  4+N:  FORMAT  (81 lo) 

FIELDS  1-8:  DEPENDENT  NODE  LABELS 

RECORDS  3,  4,  AND  4+N  ARE  REPEATED  FOR  EACH  COUPLED  DOF  OF 

THE  CONSTRAINT  SET. 

FOR  CONSTRAINT  TYPE  «  3  -  MULTI-POINT  CONSTRAINT  EQUATION 

RECORD  3:  FORMAT  (41'0, IPEli  5) 

FIELD  1:  EQUAiiC'N  .ABEL 

FIELD  2:  NUMBER  OF  TERMS 

FIELD  3:  FORCE/DISPLACEMENT  SWITCH 

1  -  FORCE 

2  -  DISPLACEMENT 


FIELD  4: 
FIELD  5: 


COLOR  NUMBER 

CONSTANT  FORCE/DISPLACEMENT  VALUE 


RECORD  4+N : 
FIELD  l: 
FIELD  2: 


FIELD  3: 


FORMAT  (110,12, 1  PE  1 3 .5) 

NODE  LABEL 

NODAL  DEGREE  OF  FREEDOM 

1  -  X 

2  -  Y 

3  -  Z 

4  -  X  ROTATIONAL 

5  -  Y  ROTATIONAL 

6  -  Z  ROTATIONAL 
CONSTRAINT  COEFFICIENT 


RECORDS  3,  4,  AND  4+N  ARE  REPEATED  FOR  EACH  MULTI -POINT 
CONSTRAINT  EQUATION  OF  THE  CONSTRAINT  SET. 

THE  ENTIRE  SET  OF  RECORDS,  INCLUDING  THE  SEPARATORS  AND  THE 
DATASET  TYPE  RECORDS,  IS  REPEATED  FOR  EACH  CONSTRAINT  TYPE  OF 
EACH  CONSTRAINT  SET  IN  THE  MODEL. 


Universal  Dataset 
Type:  89 

Name:  Restraint  set 
Status:  Current 

Written  by:  1-DEAS  Level  3  Supertab 
Read  by:  I-DEAS  Level  3  Supertab 
Revision  date:  15-JAN-i986 


RECORD  J:  FORMAT  (2110) 

FIELD  1:  RESTRAINT  SET  NUMBER 

FIELD  2:  RESTRAINT  TYPE 

1  -  NODAL  DISPLACEMENT 

2  -  NODAL  TEMPERATURE 

RECORD  2:  FORMAT  (20A2) 

FIELD  1:  RESTRAINT  SET  NAME 

FOR  RESTRAINT  TYPE  «  1  -  NODAL  DISPLACEMENT 

RECORD  3:  FORMAT  (2110,612) 

FIELD  1 :  RESTRAINT  LABEL 

FIELD  2:  COLOR  NUMBER 

FIELDS  3-8:  SWITCHES  FOR  DOFS  1-6 

0  -  OFF 
1  -  ON 

RECORD  4:  FORMAT  (1P6E13.5) 

FIELDS  1-6:  DISPLACEMENT  DOFS  1-6 

RECORDS  3  AND  4  ARE  REPEATED  FOR  EACH  NODE  IN  THE  RESTRAINT 

SET. 

FOR  RESTRAINT  TYPE  *  2  -  NODAL  TEMPERATURE 

RECORD  3:  FORMAT  (2110, IPEI3.5) 

FIELD  1:  RESTRAINT  LABEL 

FIELD  2:  COLOR  NUMBER 

FIELD  3:  TEMPERATURE  VALUE 

RECORD  3  IS  REPEATED  FOR  EACH  NODE  IN  THE  RESTRAINT  SET. 

THE  ENTIRE  SET  OF  RECORDS,  INCLUDING  THE  SEPARATORS  AND  THE 

DATASET  TYPE  RECORDS,  IS  REPEATED  FOR  EACH  RESTRAINT  TYPE  OF 

EACH  RESTRAINT  SET  IN  THE  MODEL. 


Universal  Dataset 

Type:  90 

Name:  Load  set 

Status:  Current 

Written  by:  I-DEAS  Level  3 

Read  by:  I-DEAS  Level  3 

Revision  date:  1 5- JAN- 1986 


Super tab 
Supertab 


RECORD  1: 
FIELD  1: 
FIELD  2: 


FORMAT  (2110) 

LOAD  SET  NUMBER 
LOAD  TYPE 

1  -  NODAL  FORCE 

2  -  NODAL  TEMPERATURE 

3  -  FINITE  ELEMENT  FACE  PRESSURE 

4  -  FINITE  ELEMENT  FACE  HEAT  FLUX 

5  -  FINITE  ELEMENT  EDGE  PRESSURE 

6  -  F1NTIE  ELEMENT  EDGE  HEAT  FLUX 

7  -  NODAL  HEAT  SOURCE 


RECORD  2: 
FIELD  1: 


FORMAT  (20A2) 

LOAD  SET  NAME 


FOR  LOAD  TYPE  -  1  NODAL  FORCES 


RECORD  3: 
FIELD  1: 
FIELD  2: 
FIELDS  3-8: 


FORMAT  (2110,612) 

NODAL  FORCE  LABEL 
COLOR  NUMBER 
SWITCHES  FOR  DOFS  1-6 
0  -  OFF 
I  -  ON 


RECORD  4: 
FIELD  1 
FIELD  2 
FIELD  3 
FIELD  4 
FIELD  5 
FIELD  6 


FORMAT  (1P6EI3.5) 

FORCE  FOR  DOF  1 
FORCE  FOR  DOF  2 
FORCE  FOR  DOF  3 
FORCE  FOR  DOF  4 
FORCE  FOR  DOF  5 
FORCE  FOR  DOF  6 


RECORDS  3  AND  4  ARE  REPEATED  FOR  EACH  NODAL  FORCE  OF 
THE  LOAD  SET. 


FOR  LOAD  TYPE  »  2 


NODAL  TEMPERATURES 


RECORD  3: 
FIELD  1: 
FIELD  2: 
FIELD  3: 


FORMAT  (2110, 1PE13 .5) 

NODAL  TEMPERATURE  LABEL 
COLOR  NUMBER 
TEMPERATURE  VALUE 


RECORD  3  IS  REPEATED  FOR  EACH  NODAL  TEMPERATURE  OF  THE 
LOAD  SET. 


FOR  LOAD  TYPE  -  3  FINITE  ELEMENT  FACE  PRESSURE 


RECORD  3: 

FORMAT  (2110,612) 

FIELD  1: 

FINITE  ELEMENT  LABEL 

• 

FIELD  2: 

COLOR  NUMBER 

FIELDS  3-8: 

SWITCHES  FOR  FACES  1-6 

0  -  OFF 

1  -  ON 

RECORD  4: 

FORMAT  (1P6E13.5) 

FIELDS  1-6: 

PRESSURE  VALUES  ON  FACES  1 

-6 

RECORDS  3  AND  4  ARE  REPEATED  FOR  EACH  FINITE  ELEMENT  FACE  i 

% 

PRESSURE  OF  THE  LOAD  SET. 

FOR  LOAD  TYPE  -  4 

FINITE  ELEMENT  FACE  HEAT  FLUX 

RECORD  3: 

FORMAT  (2110,612) 

FIELD  1: 

FINITE  ELEMENT  LABEL 

FIELD  2: 

COLOR  NUMBER 

FIELDS  3-8: 

SWITCHES  FOR  FACES  1-6 

0  -  OFF 

1  -  ON 

RECORD  4 : 

FORMAT  (1P6E13.5) 

FIELDS  1-6: 

HEAT  FLUX  VALUES  ON  FACES 

1-6 

RECORDS  3  AND  4  ARE  REPEATED  FOR  EACH  FINITE  ELEMENT  FACE  1 

HEAT  FLUX  OF  THE 

LOAD  SET. 

• 

FOR  LOAD  TYPE  «  5 

FINITE  ELEMENT  EDGE  PRESSURE 

RECORD  3: 

FORMAT  (2110) 

FIELD  1: 

FINITE  ELEMENT  LABEL 

FIELD  2: 

COLOR  NUMBER 

RECORD  4: 

FORMAT  (1P4E13.5) 

FIELD  1: 

PRESSURE  VALUE  ON  EDGE  1 , 

END  1 

FIELD  2: 

PRESSURE  VALUE  ON  EDGE  1 , 

END  2 

FIELD  3: 

PRESSURE  VALUE  ON  EDGE  2, 

END  1 

FIELD  4: 

PRESSURE  VALUE  ON  EDGE  2, 

END  2 

RECORD  5: 

FORMAT  (1P4E13.5) 

FIELD  1: 

PRESSURE  VALUE  ON  EDGE  3, 

END  1 

FIELD  2: 

PRESSURE  VALUE  ON  EDGE  3, 

END  2 

FIELD  3: 

PRESSURE  VALUE  ON  EDGE  4, 

END  1 

• 

FIELD  4: 

PRESSURE  VALUE  ON  EDGE  4, 

END  2 

RECORD  6: 

FORMAT  (1P4E13.5) 

FIELD  1: 

PRESSURE  VALUE  ON  EDGE  5, 

END  1 

« 

FIELD  2: 

PRESSURE  VALUE  ON  EDGE  5, 

END  2 

FIELD  3: 

PRESSURE  VALUE  ON  EDGE  6, 

END  1 

FIELD  4: 

PRESSURE  VALUE  ON  EDGE  6, 

END  2 

• 

RECORD  7: 

FORMAT  (1P4E13.5) 

FIELD  1: 

PRESSURE  VALUE  ON  EDGE  7, 

END  1 

FIELD  2 

PRESSURE  VALUE 

ON 

EDGE 

7,  END  2 

FIELD  3 

PRESSURE  VALUE 

ON 

EDGE 

8,  END  1 

FIELD  4 

PRESSURE  VALUE 

ON 

EDGE 

8,  END  2 

RECORD  8: 

FORMAT  (1P4E13.5) 

FIELD  1 

PRESSURE  VALUE 

ON 

EDGE 

9,  END  1 

FIELD  2 

PRESSURE  VALUE 

ON 

EDGE 

9,  END  2 

FIELD.  3 

PRESSURE  VALUE 

ON 

EDGE 

10,  END 

1 

FIELD  4 

PRESSURE  VALUE 

ON 

EDGE 

10,  END 

2 

RECORD  9: 

FORMAT  (1P4E13.5) 

FIELD  1 

PRESSURE  VALUE 

ON 

EDGE 

11,  END 

1 

FIELD  2 

PRESSURE  VALUE 

ON 

EDGE 

11,  END 

2 

FIELD  3 

PRESSURE  VALUE 

ON 

EDGE 

12,  END 

1 

FIELD  4 

PRESSURE  VALUE 

ON 

EDGE 

12,  END 

2 

RECORD  10 

FORMAT  (1P4E13.5) 

FIELD  1 

PRESSURE  VALUE 

ON 

EDGE 

13,  END 

I 

FIELD  2 

PRESSURE  VALUE 

ON 

EDGE 

13,  END 

2 

FIELD  3 

PRESSURE  VALUE 

ON 

EDGE 

14,  END 

1 

FIELD  4 

PRESSURE  VALUE 

ON 

EDGE 

14,  END 

2 

RECORD  11 

• 

• 

FORMAT  (1P4E13.5) 

FIELD  1 

• 

• 

PRESSURE  VALUE 

ON 

EDGE 

15,  END 

1 

FIELD  2: 

PRESSURE  VALUE 

ON 

EDGE 

15,  END 

2 

FIELD  3 

• 

PRESSURE  VALUE  ON 

EDGE 

16,  END 

1 

FIELD  4 

: 

PRESSURE  VALUE 

ON 

EDGE 

16,  END 

2 

RECORDS  3  -  1 1  ARE  REPEATED  FOR  EACH  FINITE  ELEMENT  EDGE 
PRESSURE  OF  THE  LOAD  SET.  SEE  DOCUMENTATION  FOR  THE  DEFINITIONS 
OF  EDGES  1-16  FOR  SPECIFIC  FINITE  ELEMENT  FAMILIES. 


FOR  LOAD  TYPE  -  6 

FINITE  ELEMENT  EDGE 

HEAT  FLUX 

RECORD  3: 

FORMAT  (2110,412) 

FIELD  1: 

FINITE  ELEMENT  LABEL 

FIELD  2: 

COLOR  NUMBER 

FIELDS  3-6: 

EDGE  SWITCHES  FOR  EDGES 

1  - 

4 

-0,  NO  FLUX 

ON 

THIS 

EDGE 

-1,  FLUX  ON 

THIS  EDGE 

RECORD  4: 

FORMAT  (1P4E13.5) 

FIELD  1: 

HEAT  FLUX  VALUE 

ON 

EDGE 

1, 

END 

1 

FIELD  2: 

HEAT  FLUX  VALUE 

ON 

EDGE 

1, 

END 

2 

FIELD  3: 

HEAT  FLUX  VALUE 

ON 

EDGE 

2, 

END 

1 

FIELD  4: 

HEAT  FLUX  VALUE 

ON 

EDGE 

2, 

FND 

2 

RECORD  5: 

FORMAT  (1P4E13.5) 

FIELD  1: 

HEAT  FLUX  VALUE 

ON 

EDGE 

3, 

END 

1 

FIELD  2: 

HEAT  FLUX  VALUE 

ON 

EDGE 

3, 

END 

2 

FIELD  3: 

HEAT  FLUX  VALUE 

ON 

EDGE 

4, 

END 

1 

FIELD  4: 

HEAT  FLUX  VALUE 

ON 

EDGE 

END 

2 

RECORDS  3,  4,  AND  5  ARE  REPEATED  FOR  EACH  FINITE  ELEMENT  EDGE  HEAT 
FLUX  OF  THE  LOAD  SET. 


FOR  LOAD  TYPE  »  7 


NODAL  HEAT  SOURCE 


RECORD  3:  FORMAT  (2110, 1PE13.5) 

FIELD  1:  NODAL  HEAT  SOURCE  LABEL 

FIELD  2:  COLOR  NUMBER 

FIELD  3:  HEAT  SOURCE  VALUE 

RECORD  3  IS  REPEATED  FOR  EACH  NODAL  HEAT  SOURCE  OF  THE 
LOAD  SET. 

THE  ENTIRE  SET  OF  RECORDS,  INCLUDING  THE  SEPARATORS  AND 
THE  DATASET  TYPE  RECORDS,  IS  REPEATED  FOR  EACH  LOAD  TYPE 
OF  EACH  LOAD  SET  IN  THE  MODEL. 


Definition  of  Edges  1  -  16  for  Finite  Element  Families 

1  ROD  EDGE  1:  AXIAL  DISTRIBUTED  FORCE 

EDGES  2  -  16:  NOT  APPLICABLE 

2  BEAM  EDGE  1:  AXIAL  DISTRIBUTED  FORCE 

EDGE  2:  Y  -  SHEAR  DISTRIBUTED  FORCE 
EDGE  3:  Z  -  SHEAR  DISTRIBUTED  FORCE 
EDGE  4:  DISTRIBUTED  TORQUE 
EDGE  5:  DISTRIBUTED  Y  -  BENDING  MOMENT 
EDGE  6:  DISTRIBUTED  Z  -  BENDING  MOMENT 
EDGES  7  -  16:  NOT  APPLICABLE 

3  PIPE  EDGE  1 :  AXIAL  DISTRIBUTED  FORCE 

EDGE  2:  Y  -  SHEAR  DISTRIBUTED  FORCE 
EDGE  3:  Z  -  SHEAR  DISTRIBUTED  FORCE 
EDGE  4:  DISTRIBUTED  TORQUE 
EDGE  5:  DISTRIBUTED  Y  -  BENDING  MOMENT 
EDGE  6:  DISTRIBUTED  Z  -  BENDING  MOMENT 
EDGES  7  -  16:  NOT  APPLICABLE 

4  PLANE  STRESS  EDGE  1 :  FIRST  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 

EDGE  2:  SECOND  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  3:  THIRD  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  4:  FOURTH  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  5:  FIRST  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  6:  SECOND  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  7:  THIRD  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  8:  FOURTH  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGES  9  -  16:  NOT  APPLICABLE 

5  PLANE  STRAIN  EDGE  1 :  FIRST  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 

EDGE  2:  SECOND  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  3:  THIRD  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  4:  FOURTH  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  5:  FIRST  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  6:  SECOND  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  7:  THIRD  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  8:  FOURTH  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 


EDGES  9-16:  NOT  APPLICABLE 


6  FLAT  PLATE  EDGE  1 :  FIRST  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 

EDGE  2:  SECOND  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 
EDGE  3:  THIRD  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 
EDGE  4:  FOURTH  EDGE  OF  OUT -OF -PLANE  SHEAR  LOAD 
EDGE  5:  FIRST  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGE  6:  SECOND  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGE  7:  THIRD  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGE  8:  FOURTH  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGES  9  -  16:  NOT  APPLICABLE 

7  MEMBRANE  EDGE  1 :  FIRST  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 

EDGE  2:  SECOND  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  3:  THIRD  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  4:  FOURTH  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  5:  FIRST  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  6:  SECOND  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  7:  THIRD  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  8:  FOURTH  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGES  9  -  16:  NOT  APPLICABLE 

8  AXISYMMETRIC  EDGE  1:  FIRST  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 

SOLID  EDGE  2:  SECOND  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 

EDGE  3:  THIRD  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  4:  FOURTH  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  5:  FIRST  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  6:  SECOND  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  7:  THIRD  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  8:  FOURTH  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGES  9  -  16:  NOT  APPLICABLE 

9  THIN  SHELL  EDGE  1 :  FIRST  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 

EDGE  2:  SECOND  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  3:  THIRD  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  4:  FOURTH  EDGE  OF  DISTRIBUTED  MEMBRANE  LOAD 
EDGE  5:  FIRST  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  6:  SECOND  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  7:  THIRD  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  8:  FOURTH  EDGE  OF  DISTRIBUTED  SHEAR  LOAD 
EDGE  9:  FIRST  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 
EDGE  10:  SECOND  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 
EDGE  11:  THIRD  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 
EDGE  12:  FOURTH  EDGE  OF  OUT-OF-PLANE  SHEAR  LOAD 
EDGE  13:  FIRST  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGE  14:  SECOND  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGE  15:  THIRD  EDGE  OF  DISTRIBUTED  BENDING  LOAD 
EDGE  16:  FOURTH  EDGE  OF  DISTRIBUTED  BENDING  LOAD 


10 

THICK  SHELL 

EDGES 

1  -  16: 

NOT  APPLICABLE 

11 

SOLID 

EDGES 

1  -  16: 

NOT  APPLICABLE 

12 

RIGID  BAR 

EDGES 

1  -  16: 

NOT  APPLICABLE 

13  SPRING 


EDGES  1  -  16:  NOT  APPLICABLE 


14 

DAMPER 

EDGES  1 

-  16:  NOT  APPLICABLE 

15 

GAP 

EDGES  1 

-  16:  NOT  APPLICABLE 

16 

LUMPED  MASS 

EDGES  1 

-  16:  NOT  APPLICABLE 

17 

AXISYMMETRIC 
THIN  SHELL 

EDGE  1: 
EDGE  2: 

DISTRIBUTED  MEMBRANE  LOAD 
DISTRIBUTED  SHEAR  LOAD 

EDGE  3:  DISTRIBUTED  BENDING  LOAD 
EDGES  4  -  16:  NOT  APPLICABLE 


